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1.  INTRODUCTION 


In  support  of  the  goal  to  provide  a  plasma  engineering  capability  to  the  spacecraft  community, 
the  objective  of  the  Plasma  Interactions  with  Spacecraft  contract  is  to  develop,  incorporate,  test, 
and  validate  new  algorithms  for  the  three-dimensional  plasma-environment  spacecraft 
interactions  computational  tool,  Nascap-2k.  These  algorithms  are  required  to  self-consistently 
compute  plasma  transport  and  to  model  electromagnetic  radiation  in  the  near  to  mid  field  from 
very  low  frequency  (VLF)  (3  kHz  to  30  kHz)  antennas.  The  plasma  flow  models  can  be  used  to 
address  various  plasma  engineering  concerns  including  surface  discharges  due  to  meteoroid 
impact,  and  spacecraft  contamination  due  to  electric  propulsion  plasma  plume  effects. 

Under  this  contract,  support  was  provided  for  the  DSX  (Demonstration  and  Science 
experiments)  program.  In  addition  to  program  support,  this  included  adding  capabilities  to 
Nascap-2k  to  enable  near  field  modeling  and  developing  an  algorithm  for  the  computation  of 
surface  currents  on  a  satellite  acting  as  an  antenna.  Another  task  was  the  development  of  a 
prototype  of  Nascap-2k  RealTime  for  use  by  the  SEEFS  (Space  Environmental  Effects  Fusion 
System)  program.  A  significant  effort  was  the  development  and  integration  of  a  new  database 
and  memory  manager  software  for  Nascap-2k  The  earlier  database  software,  while  more  than 
sufficient  in  its  time,  created  severe  limitations  on  further  Nascap-2k  development  and 
calculation  ability.  The  new  database  and  memory  manager  successfully  removes  these 
limitations. 

1.1  Related  Documents 

The  following  three  documents  were  prepared  under  this  contract  and  are  being  delivered  under 
separate  cover. 

N2kDB  Database  and  Memory  Manager  Software  for  Nascap-2k,  V.A.  Davis,  N.R.  Baker,  B.G. 
Gardner,  R.A.  Kuharski,  M.J.  Mandell,  A.J.  Ward,  K.G.  Wilcox 

Nascap-2k  Scientific  Documentation,  M.J.  Mandell  and  V.A.  Davis 

Nascap-2k Programmer’s  Documentation,  V.A.  Davis,  M.J.  Mandell,  K.G.  Wilcox 

1.2  Nascap-2k  Software  Development 

Nascap-2k  was  developed  by  Science  Applications  International  Corporation  (SAIC)  with 
funding  from  the  National  Aeronautics  and  Space  Administration  (NASA)  and  the  Air  Force 
Research  Laboratory  (AFRL).  Under  this  contract,  Nascap-2k  was  improved  with  new  and 
enhanced  algorithms,  and  its  calculation  capabilities  were  expanded.  Here,  we  list  many  of  the 
specific  tasks  performed  under  this  effort. 

The  most  significant  development  effort  has  been  the  design,  creation,  and  implementation  of 
N2kDB,  the  new  database  and  memory  management  software,  summarized  in  Section  1.3  below. 
This  software  is  documented  in  N2kDB  Database  and  Memory  Manager  Software  for 
Nascap-2k,  delivered  under  separate  cover. 
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A  number  of  enhancements  were  made  specifically  to  enable  DSX  near-field  calculations.  These 
are  summarized  in  Section  1 .2  below. 

At  contract  end,  in  February  201 1,  we  are  delivering  Nascap-2k  4.1  RC  (release  candidate).  This 
version  includes  the  new  database  and  memory  management  software  and  the  latest  code 
improvements.  In  January  2009,  before  incorporation  of  the  new  database  and  memory  manager 
software,  Nascap-2k  3.2  was  released  to  AFRL  and  the  NASA  SEE  Program,  which  handles 
distribution.  Nascap-2k  3.1  was  released  to  AFRL  and  the  NASA  SEE  Program  in  July  2005.  All 
of  these  versions  of  Nascap-2k  were  tested  on  our  standard  suite  of  test  cases  on  both  Windows 
and  Linux.  In  addition,  a  number  of  interim  deliveries  of  the  code  to  AFRL  were  made  (April 
2006,  July  2006,  March  2008,  October  2008,  August  2010,  and  November  2010). 

To  develop  and  test  the  parallel  processing  capabilities  and  64-bit  compatibility  of  Nascap-2k, 
we  ported  the  code  to  a  four  processor  Apple  Macintosh  Pro  with  the  MacOS  X.5  operating 
system  with  Intel  C++  and  Fortran  compilers.  To  resolve  LINUX  installation  issues,  we 
purchased  a  computer  with  a  64-bit  AMD  multicore  microprocessor  and  installed  CentOS  v5.4 
and  Portland  Group  compiler  suite  10.9.  With  this  new  computer  we  resolved  the  final  issues  and 
Nascap-2k  is  now  fully  compatible  with  the  Windows,  CentOS  LINUX,  and  64-bit  MacOS  X 
environments. 

N2kScriptRunner,  a  C++  code  that  runs  a  Nascap-2k  script  outside  of  the  Java  user  interface, 
was  created.  Using  N2kScriptRunner,  we  compiled,  linked,  and  ran  Nascap-2k  using  the 
OpenMP  compiler  commands  for  multiprocessor  operations. 

An  improvement  was  made  to  the  algorithm  used  to  display  fields  (such  as  charge  density)  that 
are  stored  as  nodal  extensive  quantities.  The  approach  used  is  described  in  Section  2. 

In  addition  to  the  above,  we  made  a  number  of  small  changes  to  Nascap-2k.  The  most  notable 
are  the  following: 

Revised  the  color  scale  used  on  the  Result  3D  tab.  The  new  color  scale  maps 
monotonically  to  gray  levels. 

Added  monopole  and  debye  screening  boundary  conditions  to  reduce  perturbations  to 
free  space  problems  caused  by  a  grounded  boundary. 

Made  a  number  of  changes  to  the  Particle  Tracking  and  Potentials  in  Space  modules 
recommended  by  AFRL  to  make  these  modules  parallelizable. 

Implemented  the  ability  to  specify  and  display  a  parameter  on  multiple  planes  at  once  on 
the  Results  3D  tab. 

Added  to  the  Results  tab  the  ability  to  view  the  surface  number  of  the  surface  with  the 
minimum  and  maximum  value  of  the  displayed  quantity. 
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Added  to  the  Results  3D  tab  the  ability  to  view  the  components  of  the  current  for  all 
timesteps. 

Added  the  ability  to  display  results  from  any  timestep  on  the  Results  3D  tab. 

Added  the  ability  to  select  a  view  direction  from  the  sun,  from  the  ram,  and  from  a  user 
specified  direction  for  the  Results  3D  tab  display. 

Added  Freja  environments  as  default  auroral  environments  on  the  Environment  tab. 

Made  the  specification  of  an  insulating  surface  as  “fixed  potential”  work  properly. 

Removed  the  use  of  XML  schema  files  for  input  validation. 

Added  a  check  to  make  sure  that  the  particle  file  filename  is  no  more  than  20  characters 
(the  most  that  can  be  read  by  the  keyword  input  routines  used  by  Tracker). 

Revised  coding  so  that  a  missing  conductor  number  is  allowed. 

Implemented  the  ability  to  close  one  project  and  open  another  project  from  the  Java  user 
interface  without  shutting  down  the  code. 

Implemented  the  buttons  on  the  Problem  tab  under  LINUX. 

Revised  the  tracking  of  the  last  directory  used.  The  file  that  contains  the  latest  directory 
used  is  now  saved  either  in  the  install  directory  or  the  user’s  home  directory, 
depending  on  the  operating  system  and  access  of  the  user  to  the  install  directory. 

We  added  a  warning  when  a  user  saves  an  object  in  which  a  node  of  one  surface  element 
is  contained  within  another  surface  element. 

The  user  interface  was  revised  to  allow  for  easy  use  of  the  new  capabilities.  The  user 
documentation  was  updated  to  reflect  the  user  interface  and  code  changes. 

We  discussed  collaborating  with  AFRL/RZSS  to  make  RZSS’s  code  COLISEUM  and  Nascap- 
2k  work  together. 

We  wrote  Nascap-2k  Scientific  Documentation  and  Nascap-2k  Programmer ’s  Documentation, 
both  delivered  under  separate  cover. 

1.3  Nascap-2k  for  DSX 

To  compute  the  sheath  structure  and  currents  about  the  DSX  VLF  antenna,  a  number  of 
capabilities  were  added  to  Nascap-2k.  This  included  improvements  to  Nascap-2k' s  surface 
charging  and  PIC  (Particle-in-cell)  computational  capabilities. 

The  ability  to  inject  macroparticles  carrying  charge  at  the  boundary  of  the  computational  space 
was  implemented.  The  model  now  has  the  ability  to  split  macroparticles  carrying  charge 
immediately  after  creation,  thus  creating  a  representation  of  the  thermal  distribution,  when  the 
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macroparticles  are  created  either  throughout  or  at  the  boundary  of  the  computational  space.  We 
also  implemented  the  ability  to  split  particles  carrying  either  charge  or  current  at  the  subgrid 
boundaries  as  needed.  A  discussion  of  these  additional  capabilities  and  our  testing  is  included  in 
Section  3. 

Upon  making  one  minor  code  change,  we  established  that  currents  computed  in  time-dependent 
Nascap-2k  calculations  can  be  used  to  compute  the  change  in  potential  of  spacecraft  surfaces. 
First,  it  was  verified  that  the  current  to  surfaces,  computed  during  tracking,  is  used  to  compute 
the  change  in  potential  and,  second,  that  the  current  is  being  used  correctly  in  the  calculation  of 
the  change  in  surface  potential.  We  also  added  the  ability  to  include  an  analytic  electron  current 
in  a  Hybrid  PIC  charging  calculation.  This  capability  is  discussed  in  Section  4. 

The  Charge  Surfaces  module  of  Nascap-2k  was  modified  so  that  the  user  may  specify  a  time- 
varying  bias  value  consisting  of  multiple  Fourier  components.  The  conductor  potentials  are 
appropriately  adjusted  to  account  for  the  internal  current  flow  as  the  bias  potential  changes. 

The  Nascap-2k  user  interface  was  modified  to  allow  the  ability  to  specify  a  loop  within  the 
script.  We  added  iteration-number-dependent  execution  of  the  Save  and  Create  Particles 
commands. 

A  self-consistent  calculation  was  performed  of  the  space  potentials  and  current  to  the  CHAWS 
(Charging  Hazards  And  Wake  Studies)  probe  on  the  Wake  Shield  Facility  (WSF)  using  the 
Hybrid  PIC  approach.  The  previous  charge  stabilization  algorithm  is  ineffective  in  dense 
plasmas.  With  the  incorporation  of  the  new  database  it  became  possible  to  save  the  volume  ion 
density  during  PIC  calculations.  This  allowed  us  to  implement  a  new  algorithm  that  linearizes 
about  the  barometric  potential  rather  than  about  zero  potential.  We  also  added  the  ability  to 
generate  a  variable  positional,  angular,  and  energy  distribution  of  macroparticles  for  PIC 
boundary  injection.  After  making  these  changes,  we  repeated  the  CHAWS  calculation  using  the 
Hybrid  PIC  approach.  The  results  are  described  in  Section  5. 

We  also  added  the  optional  ability,  when  tracking  macroparticles  carrying  charge,  to  deposit 
charge  on  the  nodes  at  the  end  of  each  sub-step  rather  than  at  the  end  of  the  time  step.  The 
applicability  and  stability  of  this  numeric  technique  for  typical  Nascap-2k  plasma  physics 
problems  is  under  evaluation.  For  a  simple  problem,  we  showed  that  as  long  as  there  are  enough 
macroparticles,  standard  PIC  and  the  orbit-averaged  approach  give  the  same  results.  This 
discussion  appears  in  Section  6. 

We  added  to  Nascap-2k  the  infrastructure  need  for  PIC  (particle-in-cell)  calculations  with  two 
species  with  different  timescales  (hydrogen  ions  and  electrons). 

A  technique  was  developed  for  the  evaluation  of  transverse  surface  currents  along  the  DSX 
antennas  and  implemented  it  in  Nas cap-2 k.  We  developed  a  related  algorithm  for  the 
computation  of  volume  electron  currents.  At  present  this  algorithm  is  implemented  in  a  stand¬ 
alone  Java  code  and  uses  ion  densities  from  an  existing  Nascap-2k  database.  A  description  of  the 
algorithms  is  in  Section  7. 
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We  added  the  ability  to  use  surface  and  volume  currents  to  compute  the  vector  potential, 
magnetic  field,  and  time  derivative  of  the  vector  potential.  The  earlier  implementation  of  volume 
ion  currents  was  corrected.  A  description  of  the  algorithm  is  in  Section  7. 

The  user  interface  was  revised  to  allow  for  easy  use  of  the  new  capabilities.  The  user 
documentation  was  updated  to  reflect  the  user  interface  and  code  changes. 

1.4  New  Database  and  Memory  Manager 

Early  in  this  contract,  we  reviewed  the  original  Nascap-2k  database  and  memory  management 
system,  the  specification,  the  desired  capabilities  of  the  replacement,  and  available  open  source 
databases  that  would  be  available  for  our  use.  We  determined  the  requirements  of  a  new  database 
and  memory  management  system  for  Nascap-2k.  The  requirements  document  appeared  as  an 
appendix  to  AFRL-VS-HA-TR-2007-1062,  the  first  interim  report  for  this  contract. 

During  the  fourth  year  of  this  contract,  we  focused  on  designing,  building,  and  implementing 
into  Nascap-2k  the  new  database  and  memory  management  software,  N2kDB.  During  the  fifth 
and  sixth  years,  we  revised  Nascap-2k  to  take  full  advantage  of  the  new  structure.  We  revised  the 
way  quantities  associated  with  conductors,  material  properties,  particle  species,  and  other  data 
are  saved  in  order  to  increase  (or  eliminate)  limitations  on  the  numbers  of  conductors,  materials, 
species,  and  timesteps. 

A  Software  Design  Document  was  produced  that  describes  the  software  design,  test  procedures, 
and  software  standards.  The  document  was  revised  during  development  as  the  design  matured. 
Versions  appeared  as  appendices  to  several  quarterly  reports.  Sections  of  this  document  were 
subsequently  reorganized  into  N2kDB  Database  and  Memory  Manager  Software  for  Nascap-2k. 

We  developed  N2kDB  Test,  a  database  testbed  code.  N2kDB  Test  has  the  same  structure  as 
Nascap-2k  and  uses  the  database  in  the  same  manner.  It  is  a  small  manageable  code  that  we  used 
for  testing  the  database  software  during  development  and  for  developing  the  proper 
implementation  of  the  new  commands.  After  incorporation  of  N2kDB  into  Nascap-2k  this  code 
is  no  longer  useful  and  is  no  longer  maintained. 

We  developed  two  versions  of  N2kDBTool  (console-based  and  with  a  Java  interface),  a  stand 
alone  program  that  reads  and  writes  Nascap-2k  database  files.  This  program  proved  invaluable 
during  development  of  N2kDB  and  promises  to  be  useful  in  future  Nascap-2k  development. 
N2kDBTool  is  included  with  Nascap-2k. 

N2kDB  also  includes  two  testing  programs  (msiotest  and  dmtest)  that  verify  the  behavior  of  the 
database  code  itself.  They  verify  that  the  requested  operation  is  performed  correctly.  (N2kDB 
Test  was  used  to  verify  that  the  appropriate  operation  is  requested.)  These  codes  will  be 
maintained  with  N2kDB,  but  will  not  be  distributed  to  users. 

N2kDB  satisfies  all  the  requirements  specified  in  the  software  requirements  document.  A  review 
of  how  these  requirements  are  satisfied  is  given  in  an  appendix  to  N2kDB  Database  and  Memory 
Manager  Software  for  Nascap-2k. 

Before  incorporation  of  N2kDB  into  Nascap-2k,  we  released  Nascap-2k  3.2. 
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In  order  to  simplify  the  transition,  we  first  implemented  N2kDB  into  Nascap-2k  using  wrappers 
between  the  original  database  commands  and  the  new  ones.  The  old  database  commands  were 
subsequently  replaced  with  direct  calls  to  the  new  database  software.  The  code  was  tested  on  a 
suite  of  problems  after  each  stage. 

1.5  Nascap-2k  RealTime 

We  developed  a  prototype  for  Nascap-2k  RealTime,  a  computer  code  that  computes  surface 
potentials  on  spacecraft  in  response  to  tabular  spectra  generated  by  magnetospheric  models.  It 
uses  a  robust  version  of  the  charging  algorithms  developed  for  Nascap-2k  and  is  an  independent 
executable  written  in  Java,  using  code  originally  developed  for  the  SEE  Spacecraft  Charging 
Handbook.  The  charging  algorithms  are  only  appropriate  to  the  plasma  environment  found  at 
geostationary  altitude  and  the  sun  direction  computation  also  assumes  that  the  spacecraft  is  at 
geostationary  altitude.  The  most  important  feature  of  this  code  is  that  it  runs  reliably  and  fast. 
Final  documentation  of  this  code  appears  in  Section  8. 

Once  Nascap-2k  RealTime  was  developed,  we  examined  the  coupling  of  a  magnetospheric 
model  ( MSM)  with  Nascap-2k  RealTime.  We  used  a  simplified  geometric  model  of  DSCS-III. 
After  several  preliminary  calculations,  we  selected  appropriate  material  properties.  Using  this 
model,  we  calculated  frame  charging  for  three  days  using  MSM  output  generated  using  three 
different  MSM  input  parameter  sets.  The  results  were  included  in  the  presentation  prepared  by 
Dr.  Robert  Hilmer  of  AFRL,  AGU  Fall  Meeting  Paper  SM41A-1 169,  “Spacecraft  Surface 
Charging  Application  Development  for  Geosynchronous  Orbit,”  R.V.  Hilmer,  D.L.  Cooke,  M. 
Tautz,  V.A.  Davis,  M.  J.  Mandell,  and  R.A.  Kuharski. 

1.6  MEO  Radiation 

We  analyzed  pitch-angle  distributions  from  the  CRRES  MEA  and  HEEF  electron  detectors.  We 
examined  the  anisotropy  factor  and  the  perpendicular  component  of  the  flux.  Details  are  given  in 
Section  9. 

1.7  DSX  Calculations 

In  2007  we  used  Nascap-2k  to  perform  a  set  of  self-consistent  calculations  of  the  plasma 
response  to  a  high  voltage  square  wave  VLF  antenna  using  the  capabilities  implemented  at  that 
time.  These  calculations  were  the  first  full  test  of  all  of  the  new  capabilities.  We  learned  that  the 
inclusion  of  a  representation  of  the  thermal  distribution  at  the  boundary  can  influence  the  current 
collected  by  the  antenna  arms.  These  calculations  are  included  in  a  paper  prepared  for  the  2007 
Spacecraft  Charging  Technology  Conference  (M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T. 
Wheelock,  C.J.  Roth,  Nascap-2k  Self-consistent  Simulations  of  a  VLF  Plasma  Antenna, 
Spacecraft  Charging  Technology  Conference,  Biarritz,  France,  June  2007). 

After  implementation  of  the  new  database  and  the  other  code  revisions  discussed  above,  we 
repeated  these  calculations.  The  new  calculations  use  an  updated  geometric  model,  have  a  larger 
computational  space,  include  a  larger  number  of  macroparticles,  and  have  improved  resolution 
about  the  antenna.  The  calculations  are  included  in  a  paper  prepared  for  the  2010  Spacecraft 
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Charging  Technology  Conference  (M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  C.J. 
Roth,  Nascap-2k  Self-consistent  Simulations  of  a  VLF  Plasma  Antenna,  Spacecraft  Charging 
Technology  Conference,  Albuquerque,  New  Mexico,  September  2010). 

Since  September  2010,  we  extended  the  calculations,  performed  calculations  including  any  effect 
of  the  magnetic  field  on  the  ion  motion,  and  performed  a  calculation  at  higher  voltage.  There  is 
no  magnetic  field  effect  on  ion  motion.  Increasing  the  applied  voltage  from  1  kV  to  5  kV  led  to 
no  apparent  qualitative  differences. 

The  final  results  are  discussed  in  Section  10. 

1.8  DSX  Program  Support 

In  addition  to  calculations  of  the  near  field  response  to  the  DSX  antenna,  we  supported  the  DSX 
program  by  participating  in  various  meetings  and  teleconferences. 

Via  teleconference,  Dr.  Mandell  attended  the  October  4,  2006  DSX  HSD  Solar  Array  Subsystem 
Preliminary  Design  Review  held  in  Littleton,  CO. 

We  took  an  active  role  in  the  Y-boom  high  voltage  isolation  design  discussion.  To  that  end,  Dr. 
Myron  J.  Mandell  participated  in  the  following  teleconferences  and  meetings. 

•  Weekly  teleconferences  with  Y  Antenna  supplier  (L’Garde)  from  November  6,  2007  to 
through  February  13,  2008. 

•  A  visit  to  the  L’Garde  facility  to  discuss  the  design  issues  on  November  7,  2007. 

•  Y  Antenna  Isolation  teleconference  on  December  19,  2007. 

•  Weekly  teleconferences  with  the  alternate  supplier  of  Y  Antenna  (ATK)  from  January  15 
to  February  25,  2008. 

•  ATK  Isolation  face-to-face  at  Houston  Airport  (IAH)  on  February  13,  2008. 

•  DSX  Alternative  Y  Antenna  CDR,  Goleta,  California  on  February  26  by  teleconference. 

•  DSX  Y  Antenna  CDR  at  L’Garde,  Tustin,  California  on  February  27. 

Dr.  Myron  Mandell  attended  and  made  a  presentation  at  the  Workshop  on  The  Remediation  of 
Enhanced  Radiation  Belts  in  Lake  Arrowhead,  California  on  March  3-6,  2008. 

Dr.  Mandell  attended  the  DSX  System  CDR,  Breckenridge,  Colorado,  May  6-8,  2008. 

Dr.  Mandell  attended  the  MURI  Review  and  RBR  Workshop  at  Stanford  in  Palo  Alto, 

California,  February  18-19,  2009. 

Dr.  Mandell  and  Dr.  Davis  attended  the  DSX  Science  Team  Meeting  at  Lake  Arrowhead,  CA,  on 
September  15-18,  2009.  Dr.  Mandell  made  a  presentation. 
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1.9  Contract 


The  scientists  and  other  researchers  who  contributed  to  this  work  are  Dr.  Myron  J.  Mandell,  Dr. 
Victoria  A.  Davis,  Dr.  Stuart  L.  Huston,  Dr.  Robert  A.  Kuharski,  Dr.  Michael  Brown-Hayes,  Ms. 
Barbara  M.  Gardner,  Ms.  Katherine  Wilcox,  Ms.  Alisa  J.  Ward,  and  Mr.  Nicholas  R.  Baker. 

This  contract  continues  work  performed  under  earlier  contracts:  F 1 9628-9 l-C-0 187,  Space 
System-Environment  Interactions  Investigation;  F19628-93-C-0050,  Modeling  and  Post  Mission 
Data  Analysis;  F19628-89-C-0032,  Analysis  of  Dynamical  Plasma  Interactions  with  High 
Voltage  Spacecraft;  and  F19628-98-C-0074,  Spacecraft  Potential  Control.  NASA  supported 
related  work  under  contracts  NAS8-98220  and  NAS8-02028. 

1.10  Publications 

The  following  publications  were  supported  in  total  or  in  part  by  this  contract. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  Nascap-2k  Spacecraft  Charging  Code 
Overview,  Proceedings  of  the  9th  Spacecraft  Charging  Technology  Conference,  Tsukuba,  Japan, 
2005. 

V.A.  Davis,  M.J.  Mandell,  F.J.  Rich,  D.L.  Cooke,  Reverse  trajectory  approach  to  computing 
ionospheric  currents  to  the  Special  Sensor  Ultraviolet  Limb  Imager  on  DMSP,  Proceedings  of 
the  9th  Spacecraft  Charging  Technology  Conference,  Tsukuba,  Japan,  2005. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Nascap-2k  simulations  of  a 
VLF  plasma  antenna,  Proceedings  of  the  9th  Spacecraft  Charging  Technology  Conference, 
Tsukuba,  Japan,  2005. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  Nascap-2k  spacecraft  charging  code 
overview,  IEEE  Trans  Plasma  Science,  34,  No.  5,  p  2084,  2006. 

V.A.  Davis,  M.J.  Mandell,  F.J.  Rich,  D.L.  Cooke,  Reverse  trajectory  approach  to  computing 
ionospheric  currents  to  the  Special  Sensor  Ultraviolet  Limb  Imager  on  DMSP,  IEEE  Trans 
Plasma  Science,  34,  No.  5,  p  2062,  2006. 

V.A.  Davis,  M.J.  Mandell,  D.L.  Cooke,  D.C.  Ferguson,  Nascap-2k  spacecraft  plasma 
environment  interactions  modeling:  Capabilities  and  verification,  AIAA  2007-1096,  Aerospace 
Sciences  Meeting  and  Exhibit,  Reno,  2007. 

M.J.  Mandell,  V.A.  Davis,  B.M.  Gardner,  F.K.  Wong,  R.C.  Adamo,  D.L.  Cooke,  A.T. 

Wheelock,  Charge  Control  of  Geosynchronous  Spacecraft  using  Field  Effect  Emitters,  AIAA 
2007-284,  Aerospace  Sciences  Meeting  and  Exhibit,  Reno,  2007. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Nascap-2k  Self-consistent 
Simulations  of  a  VLF  Plasma  Antenna,  Spacecraft  Charging  Technology  Conference,  Biarritz, 
France,  June  2007. 
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M.J.  Mandell,  V.A.  Davis,  E.J.  Pencil,  M.J.  Patterson,  H.K.  McEwen,  J.E.  Foster,  J.S  Snyder, 
Modeling  the  NEXT  Multi-Thruster  Array  Test  with  Nascap-2k,  Spacecraft  Charging 
Technology  Conference,  Biarritz,  France,  June  2007.  (Original  research  supported  by  NASA.) 

M.J.  Mandell,  V.A.  Davis,  E.J.  Pencil,  M.J.  Patterson,  H.K.  McEwen,  J.E.  Foster,  J.S  Snyder, 
Modeling  the  NEXT  Multi-Thruster  Array  Test  with  Nascap-2k,  IEEE  Transactions  on  Plasma 
Science,  36,  p.  2309,  2008.  (Original  research  supported  by  NASA.)  V.A.  Davis,  M.J.  Mandell, 
D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Nascap-2k  Self-consistent  Simulations  of  a  VLF  Plasma 
Antenna,  Spacecraft  Charging  Technology  Conference,  Albuquerque,  New  Mexico,  September 
2010. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Pseudopotential  algorithms 
for  simulation  of  VLF  plasma  antenna  current  flow,  Spacecraft  Charging  Technology 
Conference,  Albuquerque,  New  Mexico,  September  2010. 

We  made  presentations  or  supported  others  who  made  presentations  at  the  following  meetings. 

9th  Spacecraft  Charging  Technology  Conference,  Tsukuba,  Japan,  April  2005. 

(Papers  listed  above.) 

Sleight  of  HAND  (RBR  Phase  2  Kickoff)  Meeting,  Hanscom  AFB,  September  2005. 

M.J.  Mandell,  Antenna  Sheath  Modeling. 

AGU  Fall  Meeting,  San  Francisco,  CA,  December  2005 

R.V.  Hilmer,  D.L.  Cooke,  M.  Tautz,  V.A.  Davis,  M.  J.  Mandell,  and  R.A.  Kuharski, 
Spacecraft  surface  charging  application  development  for  geosynchronous  orbit,  SM41A- 
1169. 

45th  Aerospace  Sciences  Meeting,  Reno,  NV,  January  2007. 

(Papers  listed  above.) 

10th  Spacecraft  Charging  Technology  Conference,  Biarritz,  France,  June  2007. 

(Papers  listed  above.) 

Workshop  on  Remediation  of  Enhanced  Radiation  Belts,  Lake  Arrowhead,  CA,  March  2008. 

M.  J.  Mandell,  V.  A.  Davis,  D.  L.  Cooke,  A.  T.  Wheelock,  C.  J.  Roth,  Nascap-2k 
simulations  of  static  and  dynamic  sheaths. 

GOMACTech-08,  Las  Vegas,  NV,  March,  2008. 

M.J.  Mandell,  V.A.  Davis,  A.  Wheelock,  D.L.  Cooke,  Modeling  Space  Weather  Effects 
using  Nascap-2k. 

47th  AIAA  Aerospace  Sciences  Meeting,  Orlando,  FL,  January  2009 

D.L.  Cooke,  A.T.  Wheelock,  V.A.  Davis,  M.J.  Mandell,  Simulation  of  auroral  charging 
in  the  DMSP  environment,  AIAA  2009-350. 
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DSX  Science  Team  Meeting,  Lake  Arrowhead,  CA,  September  2009 

M.  J.  Mandell  and  V.  A.  Davis,  Nascap-2k  surface  current  calculations. 

1 1th  Spacecraft  Charging  Technology  Conference,  Albuquerque,  NM,  September  2010. 

(Papers  listed  above.) 

2.  PLOTTING  AND  POTENTIAL  CALCULATION  IMPROVEMENTS 
2.1  Monopole  Boundary  Condition 

When  potentials  are  long-range,  the  shape  and  size  of  the  sheath  may  be  substantially  affected  by 
a  zero  potential  outer  boundary.  This  was  recognized  in  NASCAP/GEO,  where  a  1/r  boundary 
potential  could  be  set  with  magnitude  proportional  to  the  computed  charge  on  the  object.  No 
such  option  was  implemented  for  subsequent  codes,  as  their  intended  regime  of  applicability  was 
dense  plasma,  and  thus  potentials  were  always  short-range. 

In  Nascap-2k  the  problem  was  often  apparent  during  tenuous  plasma  calculations.  Addition  of 
outer  grids  was  not  always  sufficient  to  adequately  minimize  the  problem,  while  adding 
substantially  to  the  problem’s  size  and  complexity. 

The  algorithm  implementing  monopole  boundaries  is  based  on  adding  finite  elements  extending 
from  the  problem  boundary  to  infinity.  The  consequent  residuals  to  be  added  to  the  boundary 
nodes  can  be  calculated  analytically  for  Laplacian  potentials.  The  algorithm  is  implemented  as 
well  for  linearly  screened  (finite  debye  length)  potentials,  although  for  that  case  some  of  the 
integrals  have  been  approximated.  While  this  treatment  makes  the  approximation  that  the  center 
of  charge  is  at  the  grid  center,  it  does  not  depend  on  knowledge  of  the  object  charge,  nor  does  it 
require  the  boundary  potential  to  be  symmetric  about  the  grid  center. 

Examples  of  potential  calculations  with  zero  potential  and  with  monopole  boundary  conditions 
are  shown  in  Figure  1. 


Figure  1.  Potential  Calculation  Using  Zero  Potential  Boundary  Conditions  (Left)  and 

Monopole  Boundary  Conditions  (Right) 
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2.2  Multiple  Cut  Planes 

We  implemented  machinery  and  a  user  interface  to  display  multiple  cut  planes  simultaneously. 
An  arbitrary  number  of  cut  planes  can  be  defined  (by  normal  direction  and  offset)  and 
individually  shown  and  hidden.  An  example  is  shown  in  Figure  2. 


Figure  2.  Example  of  Multiple  Cut  Plane  Display 


2.3  Nodal  Charge  Density 

Using  an  earlier  version  of  Nas cap-2 k,  attempts  were  made  to  plot  the  value  of  the  Nodal  Charge 
Density  on  a  cutplane.  The  plots  had  non-physical  structures,  particularly  near  grid  boundaries. 
Also,  the  Nodal  Charge  Density  on  a  plane  as  displayed  on  the  3-D  Results  tab  showed  non¬ 
physical  values  along  grid  boundaries.  The  reason  for  the  non-physical  values  is  that  the  Nodal 
Charge  Density  as  displayed  is  actually  the  charge  on  each  node  as  stored  in  the  database  and 
used  in  the  calculations  (an  extensive  quantity)  converted  to  the  charge  density  (an  intensive 
quantity)  for  the  purpose  of  display.  Extensive  quantities  exist  only  on  nodes,  and  therefore  it 
does  not  make  sense  to  discuss  their  values  at  arbitrary  points.  The  correct  way  to  convert  an 
extensive  quantity  to  an  intensive  quantity  is  essentially  to  multiply  by  the  inverse  of  the  volume 
matrix.  Unfortunately,  while  straightforward,  implementing  this  is  quite  complex.  With  a 
moderate  level  of  effort  we  were  able  to  implement  an  improvement  to  the  display  of  nodal 
extensive  quantities. 
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We  added  two  grid  interface  operations  to  the  cut  plane  display  routines.  One  grid  interface 
operation  is  before  and  one  after  the  division  of  each  node  by  mesh  volume.  In  the  first  grid 
interface  operation,  the  extensive  quantity  on  nodes  on  grid  boundaries  is  transferred  from  the 
inner  (finer)  grid  to  the  outer  (coarser)  grid.  In  the  second,  the  intensive  quantity  (resulting  from 
the  division  by  volume)  is  interpolated  from  the  outer  grid  to  the  inner  grid  boundary  nodes.  No 
adjustment  is  made  to  account  for  the  fact  that  the  appropriate  volume  for  grid  interface  nodes 
and  those  near  the  object  is  not  the  cube  of  the  mesh  unit. 

The  nodal  charge  density,  derived  from  the  nodal  charge  using  the  old  and  new  approaches  for 
display,  is  shown  in  the  Figure  3  and  Figure  4  for  the  moving  sphere  problem.  Most  of  the 
artificial  structures  on  the  grid  boundaries  have  been  eliminated. 


Figure  3.  Nodal  Charge  Density  on  X  =  0  Plane  Using  Original  and  Revised  Algorithm  for 

Display  of  Extensive  Quantity 


12 


Density 

Volts/m2 
26.0 


Before 


After 


Density 

Volts/m2 
26.0 


Figure  4.  Nodal  Charge  Density  on  Z  =  0  Plane  Using  Original  and  Revised  Algorithm  for 

Display  of  Extensive  Quantity 

3.  SPLITTING  AND  INJECTING  PARTICLES  IN  NASCAP-2K 

This  section  outlines  coding  and  testing  performed  for  splitting  of  particles  in  Nascap-2k.  The 
desirability  of  particle  splitting  has  become  apparent  both  to  avoid  having  heavy  particles  in 
well-resolved  regions  and  to  simulate  a  thermal  distribution.  Splitting  is  done  in  such  a  way  that 
merging  of  particles  is  straightforward,  although  it  is  not  obvious  that  particle  merging  is  needed 
in  Nascap-2k.  Injecting  thermal  particles  at  boundaries  is  also  part  of  this  effort. 

3.1  General  Principles 

1 .  Particles  are  split  in  velocity  space  only.  Because  we  frequently  find  ourselves  in  high- field 
regions,  spatial  splitting  would  raise  problems  with  energy  conservation. 

2.  To  be  split  in  velocity  space,  each  particle  must  carry  a  temperature.  We  assume  the 
temperature  is  always  isotropic.  The  fission  products  carry  half  the  temperature  of  the 
original  particle,  while  the  remaining  thermal  energy  appears  as  kinetic  energy  of  the  split 
particles.  PartGenDLL  initializes  the  temperature  to  the  value  of  TION  (even  for  electrons). 

3.  For  splitting  purposes,  we  define  the  Z-axis  to  be  along  the  direction  of  the  particle  velocity, 
the  X-axis  randomly  chosen  in  the  plane  normal  to  Z,  and  the  Y-axis  mutually  perpendicular. 

4.  We  split  into  two  or  three  particles  with  added  velocity  along  each  axis,  except  that  we  may 
elect  not  to  split  along  the  Z-direction  if  the  kinetic  energy  exceeds  the  thermal  energy.  Not 
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splitting  along  Z  helps  ameliorate  particle  proliferation,  but  makes  an  error  by  not  preserving 
the  original  particle  temperature  along  Z.  We  thus  end  up  with  eight,  nine,  or  twenty-seven 
new  particles. 

5.  Particle  velocity  is  assumed  to  be  acquired  by  acceleration  rather  than  actual  drift  (i.e., 
spacecraft  velocity).  If  there  is  actual  drift  (e.g.,  ram  velocity),  then  the  drift  velocity  should 
be  removed  before  splitting  the  particle,  and  added  back  after. 


6. 


If  the  particle  with  speed  uo  is  split  by  two  along  the  X  or  Y  axis,  the  new  velocity  is 
±0.707^T/m  Along  the  Z  axis,  the  velocity  increment  is  calculated  as  if  the  temperature 
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7.  If  the  particle  with  speed  uo  is  split  by  three  along  the  X  or  Y  axis,  there  is  a  zero-velocity 
central  particle  and  two  “probe”  particles  with  velocity  ±0.S66^T/m  .  Along  the  Z  axis,  the 


velocity  increment  is  calculated  as  if  the  temperature  were  T 
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3.2  Implementation 

Particle  splitting  has  been  implemented  in  PartGenDLL  for  particles  read  from  an  external  file, 
for  space-filling  default  particles,  and  for  particles  injected  from  the  boundary.  Particles  are  split 
if  the  keyword  SPLIT  appears  in  the  input  file.  Particles  read  from  an  external  file  are  split  using 
the  default  option  and  the  others  are  split  using  the  eight-particle  option. 

3.3  Examples 

This  example  was  run  using  the  “Sphere  Hybrid  PIC”  problem  that  we  have  been  using  to 
develop,  document,  and  test  PIC  currents  and  charging.  A  2.4-m  cubic  space  is  filled  with  H+ 
with  1-eV  temperature  and  density  10  m'  .  The  ions  are  collected  by  a  10-cm  radius  probe 
biased  to  -100  V.  The  collection  of  ions  by  the  probe  and  the  loss  of  ions  out  the  sides  are 
monitored,  and  the  final  potential  and  particle  configurations  are  inspected.  Conceptual  errors 
were  found  with  the  way  this  problem  was  set  up,  mostly  relating  to  the  velocity  initialization  by 
subroutine  INIVEL  prior  to  the  splitting  of  the  particles. 

Figure  5  shows  the  initial  results,  with  the  calculation  performed  with  the  old  velocity 
initialization,  with  the  addition  of  particle  splitting.  The  collected  current,  estimated  to  be  40  pA 
in  equilibrium,  rapidly  rises  to  a  sustained  value  of  about  100  pA.  The  escaping  current, 
estimated  at  200  pA,  averages  to  a  mere  15  pA.  Figure  6  shows  the  potentials  after  25  ps.  The 
admirably  spherical  sheath  is  surrounded  by  a  ring  of  positive  potentials  (~0.3  V). 
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Figure  5.  Collected  Current  (Left  Scale)  and  Escaping  Current  (Right  Scale)  Using  Default 
Script  and  Original  INIVEL  Velocity  Initialization 


Figure  6.  Potentials  at  25  ps  Corresponding  to  Figure  5 
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The  reason  for  this  behavior  is  that  particles  in  the  field- free  region  (which  does  have  some  small 
inward  field)  are  initialized  with  inward  velocities  comparable  to  the  thermal  speed. 

Convergence  of  particles  moving  inward  through  the  field-free  region  causes  the  ring  of  positive 
potentials,  and  also  explains  the  higher  potentials  toward  the  comers.  Also,  since  particles  are 
moving  inward,  there  is  little  escaping  current. 

We  changed  the  code  to  assign  small  initial  velocities  in  regions  where  the  fields  are  small  but 
non-vanishing,  while  making  no  change  at  or  within  the  sheath.  The  results  are  shown  in  Figure 
7  and  Figure  8.  The  initial  current  rise  remains  as  before,  because  ions  within  the  sheath  were 
initialized  in  the  same  way.  However,  the  current  drops  to  a  value  of  about  60  pA,  which  is  much 
closer  to  the  analytic  estimate  of  40  pA,  especially  if  we  allow  for  presheath  enhancement.  The 
escaping  current  averages  to  about  150  pA,  which  is  acceptably  close  to  the  analytic  estimate  of 
200  pA.  The  positive  potential  region  is  gone  because  we  have  not  assigned  convergent 
velocities  to  the  particles  in  the  field-free  region.  Instead,  the  initially  field- free  region  has 
attained  a  negative  potential  of  about  -0.3  V  as  the  ion  population  is  depleted  both  by  being 
collected  and  by  escaping  the  grid. 


Time 


Figure  7.  Collected  Current  (Left  Scale)  and  Escaping  (Lost)  Current  (Right  Scale)  Using 
Default  Script  and  Modified  INIVEL  Velocity  Initialization 
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Figure  8.  Potentials  at  25  (is  Corresponding  to  Figure  7 

Finally,  in  order  to  perform  the  calculation  intended,  we  modified  the  start  of  the  calculation  to 
the  following: 

1 .  Initialize  the  probe  to  zero  potential. 

2.  Calculate  potentials  and  fields  (immediately  converging  to  all  zeros). 

3.  Create  (and  split)  particles.  The  created  particles  now  have  only  the  velocities  that  result 
from  the  splitting. 

4.  Reset  the  probe  to  -100  V. 

5.  Recalculate  potentials  (in  Hybrid  PIC  mode  using  the  created  particles). 

6.  Start  tracking. 

Results  from  this  calculation  are  shown  in  Figure  9  and  Figure  10.  The  results  differ  from  those 
of  the  previous  calculation  only  in  the  slower  rise  time  and  lower  peak  in  the  current,  which 
occurs  because  now  the  ions  must  be  accelerated  before  they  can  be  collected.  Also,  the  collected 
current  is  less  noisy. 
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Figure  9.  Collected  Current  (Left  Scale)  and  Escaping  Current  (Right  Scale)  after 
Initializing  Velocities  by  Thermal  Splitting  Only 


Figure  10.  Potentials  at  25  ps  Corresponding  to  Figure  9 
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3.4  Boundary  Injection 

In  Figure  1 1,  we  note  that  the  field-free  region  has  negative  potential  mostly  in  the  range  of  0.1 
to  0.3  V.  Presumably,  continued  collection  and  escape  of  ions  will  lower  this  potential  further. 
Boundary  injection  should  keep  these  potentials  near  zero  by  replenishing  the  ions  that  have  been 
collected  or  escaped.  We  implement  boundary  injection  with  a  new  “Create  Particles”  call  at 
each  timestep  following  the  potential  solution  and  preceding  the  particle  tracking.  The  particle 
type  is  “INJECT”  and  the  time  interval  corresponding  to  the  injection  (equal  to  the  timestep  if 
injection  is  done  every  timestep)  is  required  in  the  third  field  of  the  “PARTTYPE”  input  line. 

As  a  side  effect,  the  “Create  Particles”  call  results  in  pruning  of  the  dead  and  escaped  particles, 
so  that  injection  does  not  necessarily  result  in  increased  particle  number. 

The  implementation  is  to  have  an  injection  point  at  each  quarter-boundary-surface-element.  No 
particle  is  injected  if  the  electric  field  directs  the  particle  back  towards  the  boundary.  Otherwise, 
the  injected  particle  has  charge  equal  to  the  plasma  thermal  current  times  the  area  times  the  time 

l2eT 

interval,  and  velocity  equal  to_ . - ,  so  that  it  represents  a  density  of  n/2.  Optionally,  the 

A1  ran 

injected  particles  can  be  split  into  eighths. 

Figure  1 1  shows  the  collected,  escaping,  and  injected  current.  The  collected  and  escaping 
currents  are  indistinguishable  from  those  seen  in  Figure  9  without  boundary  injection; 
undoubtedly,  there  would  be  some  divergence  if  the  problem  was  run  longer.  The  injected 
current  is,  on  average,  slightly  greater  than  the  escaping  current,  and  far  less  than  the  sum  of  the 
collected  and  escaping  currents. 


Figure  11.  Collected  Current  (Left  Scale),  Escaping  Current  (Right  Scale),  and  Injected 
Current  (Right  Scale)  Running  Problem  with  Boundary  Injection 
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Figure  12  shows  the  potentials  in  the  presence  of  boundary  injection.  The  outermost  “clean” 
spherical  contour  is  at  -0.3  Y,  the  same  as  in  Figure  10.  In  addition,  we  have  a  somewhat  ragged 
spherical  contour  at  -0.1  V,  and  a  very  ragged  (but  still  recognizable)  contour  at  -0.03  V.  The  noise 
in  the  field-free  region  is  well  under  ±0. 1  V,  and  mostly  under  ±0.03  V.  By  contrast,  without 
boundary  injection  (see  Figure  10)  the  potential  was  more  negative  than  -0.1  V  in  nearly  all  the 
field-free  region,  with  islands  more  negative  than  -0.3  V  near  the  comers. 


Figure  12.  Potentials  at  25  ps  Corresponding  to  Figure  9  (Compare  with  Figure  8) 

We  next  experimented  with  splitting  the  boundary  injected  particles.  Figure  13  compares  the 
results  using  split  and  unsplit  injected  ions,  ran  for  50  ps.  The  collected  current  shows  no 
significant  difference.  However,  the  escaping  current  at  late  times  is  significantly  higher  (and 
thus  closer  to  the  expected  value  of  200  pA)  when  the  particles  are  split. 

The  difference  between  split  and  unsplit  injection  is  more  apparent  in  the  final  potential 
contours.  Split  injection  (Figure  15)  shows  a  transition  from  a  spherical  contour  at  -0.3  V  to  a 
square  contour  at  -0.1  V,  and  thence  a  nearly  noiseless  path  to  the  problem  boundary.  Without 
splitting  (Figure  14),  the  -0.3  volt  contour  is  already  beginning  to  square,  and  from  there  to  the 
boundary,  there  is  very  noticeable,  albeit  low-level,  noise. 

The  minor  improvement  resulting  from  particle  splitting  comes  at  considerable  cost.  From  an 
initial  particle  count  of  2.4  million,  the  ran  with  unsplit  boundary  injection  has  its  particle  count 
decrease  to  1.5  million  at  50  ps,  while  the  ran  with  split  boundary  injection  sees  an  increase  to 
8.1  million  at  50  ps. 
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Figure  13.  Collected  (left  scale)  and  Escaping  (Right  Scale)  Currents  for  Calculations  in 
which  the  Boundary  Injected  Ions  are  Split  or  Unsplit 
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Figure  15.  Potential  Contours  after  50  ps  for  Split  Boundary  Injected  Ions 
3.5  Splitting  on  Entering  a  Finer  Grid 

The  real  reason  for  assigning  temperatures  to  particles  is  so  that  they  can  be  split  repeatedly  in 
mid-flight.  Figure  16  shows  a  quadrant  of  particles  after  9  ps.  The  particles  are  initially  unsplit, 
and  the  code  parameters  (see  below)  are  set  such  that  ions  will  be  split  into  nine  or  twenty-seven 
particles  on  entering  a  finer  grid.  In  Figure  16,  ions  can  be  seen  entering  Grid  2  from  Grid  1  at 
the  top  and  right.  Because  these  ions  are  moving  slowly,  they  are  split  both  along  and  normal  to 
their  motion  direction.  By  this  time,  all  ions  originating  in  Grid  3  have  been  “eaten”  by  the 
sphere,  so  that  the  cloud  of  ions  currently  in  Grid  3  has  entered  from  Grid  2  and  been  split.  Those 
that  entered  most  recently  were  already  drifting  significantly  when  they  were  split,  and  were  thus 
split  only  normal  to  their  direction  of  motion,  as  described  above. 
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Figure  16.  Particle  Positions  after  9  jis  when  Particles  are  Split  on  Entering  a  More  Finely 

Resolved  Grid 


3.6  Moving  Frames 

It  is  important  to  be  able  to  split  and  inject  particles  for  cases  when  the  spacecraft  is  moving 
through  the  plasma.  The  calculation  is  performed  in  the  spacecraft  frame  (plasma  moving  with 
“ram  flow”),  while  the  temperature  is  measured  in  the  plasma  frame.  Therefore,  to  split  particles 
requires  transforming  the  velocity  from  the  spacecraft  frame  to  the  plasma  frame,  applying  the 
splitting  algorithm  outlined  above,  and  re-transforming  velocities  of  split  particles  back  to  the 
plasma  frame. 

When  injecting  particles,  we  compare  the  inward  component  of  the  ram  velocity  with  the  usual 

|2eT 

injection  speed,  i -  .  If  the  ram  component  is  greater  than  the  usual  injection  speed,  then  we 

\  Kill 

always  inject  with  the  ram  velocity.  Otherwise,  the  injection  velocity  is  determined  by  adding 
the  stationary  injection  velocity  to  the  ram  velocity;  if  this  is  inward,  then  we  inject  only  when 
the  electric  field  is  attractive.  The  weight  of  the  injected  particles  is  calculated  from  the  inward 
normal  component  of  the  vector  sum  of  the  ram  current  and  the  thermal  current.  The  current  and 
velocity  are  related  in  such  a  way  that  the  density  contribution  of  the  injected  particles  varies  from 
half  the  ion  density  for  a  stationary  object  to  the  full  ion  density  for  a  high  mach  number  object. 
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Figure  17  shows  a  hybrid  PIC  calculation  (with  no  particle  splitting)  of  0+  flow  past  an 
uncharged  sphere  moving  at  LEO  velocity  in  the  (1,1,0)  direction.  (Geometry  and  other 
parameters  are  the  same  as  used  thus  far  in  the  document.)  Clearly  seen  are  the  negative 
potentials  in  the  object  wake  and  the  boundary  between  the  injected  particles  (diagonal  pattern) 
and  the  original  particles  (square  pattern).  Potential  fluctuations  on  the  order  of  0.05  V  are  seen 
in  the  first  subdivided  grid  where  it  is  populated  by  outer  grid  ions. 


Figure  17.  Potentials  and  Ion  (0+)  Macroparticles  after  80  ps  for  an  Uncharged  Sphere 
Moving  in  the  (1,1?0)  Direction  (No  Splitting  of  Macroparticles) 


Figure  18  shows  a  similar  calculation,  now  with  the  sphere  once  again  charged  to  -100  V.  After 
80  ps,  ions  are  focused  in  the  wake  with  sufficient  strength  to  create  positive  potentials  of 
approximately  one  volt.  This  structure  persists,  as  shown  in  Figure  19,  where  the  potential 
maximum  in  the  wake  has  reached  nearly  the  ram  energy. 
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Figure  18.  Potentials  and  Ion  (0+)  Macroparticles  after  80  ps  for  a  Sphere  Charged  to 
-100  V  Moving  in  the  (1,1,0)  Direction  (No  Splitting  of  Macroparticles) 


Figure  19.  Same  Calculation  as  Figure  18  after  136  ps 
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Figure  20  shows  the  charged  sphere  calculation,  now  with  particles  split  on  entering  a  refined 
grid.  While  the  general  character  of  the  result  is  the  same,  the  potentials  are  much  smoother  and 
now  show  compression  of  the  sheath  on  the  ram  side.  Figure  21  shows  the  current  to  the  sphere. 
After  an  early  peak  to  nearly  16  pA,  the  current  settles  down  to  a  value  of  fewer  than  7  pA, 
comfortably  less  than  the  orbit-limited  value  of  8  pA.  Of  course,  this  improved  fidelity  comes  at 
a  price,  with  a  final  particle  count  in  excess  of  two  million  in  Figure  20,  versus  under  0.4  million 
in  Figure  19. 


Figure  20.  Potentials  and  Ion  (0+)  Macroparticles  after  160  ps  for  a  Sphere  Charged  to 
-100  V  Moving  in  the  (1,1,0)  Direction  -  Particles  Split  on  Entering  Refined  Grid 
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Time 


Figure  21.  Current  for  Sphere  Charged  to  -100  V  Moving  in  the  (1,1,0)  Direction  - 

Particles  split  on  Entering  Refined  Grid 

4.  PARTICLE-IN-CELL/CHARGING  CALCULATIONS  IN  NASCAP-2K 

We  made  several  modifications  to  Nascap-2k  to  simplify  and  enhance  time-dependent 
calculations  in  which  the  tracked  current  is  used  to  compute  the  change  in  potential  of  spacecraft 
surfaces. 

New  options  have  been  added  to  the  Problem  tab.  Time-dependent  problems  may  be  specified 
as  either  “Fixed  Surface  Potentials”  or  “Self-Consistent  Surface  Potentials.”  If  “Self-Consistent 
Surface  Potentials”  is  chosen,  either  “Tracked  Particle  Currents  Only”  or  “Tracked  Ion  & 
Analytic  Electron  Currents”  can  be  selected.  “Time-Dependent  Plasma”  calculations  with  “Fixed 
Surface  Potentials”  are  calculations  of  surface  currents  as  a  function  of  time  for  user  specified 
surface  potentials.  “Time-Dependent  Plasma”  calculations  with  “Self-Consistent  Surface 
Potentials”  are  the  same,  with  the  addition  of  a  charging  step.  The  surface  charging  is  computed 
either  with  only  the  tracked  current  (“Tracked  Particle  Currents  Only”)  or  with  both  the  tracked 
current  and  an  analytic  electron  current  (“Tracked  Ion  &  Analytic  Electron  Currents”).  If  the 
second  option  is  chosen,  then  an  environment  is  set  in  the  charging  calculation  and  an  electron 
current  that  is  a  function  of  the  surface  area,  the  electron  thermal  current,  the  surface  potential, 
and  the  plasma  temperature  is  added  to  each  surface  when  computing  the  potential. 


1  = 


Ajtli  exp(<f>/0) 
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i+^ 

0 


if  <|)<0 
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[1] 


27 


If  “Self-Consistent  Surface  Potentials”  is  chosen,  the  default  script  is  as  follows  and  is  similar  to 
that  used  for  a  LEO  charging  problem.  The  SetEnvironment  command  is  only  included  for 
“Tracked  Ion  &  Analytic  Electron  Currents.”  The  ZeroCurDeriv Algorithm  command  specifies 
that  the  change  in  potential  due  to  the  tracked  current  in  a  single  timestep  is  calculated  using  an 
explicit  algorithm  rather  than  an  implicit  one.  The  change  in  the  current  due  to  the  change  in 
potential  during  the  timestep  is  not  included  when  calculating  the  potential  change.  The 
DoTrackTimeStep  command  is  the  same  as  the  DoOneTimeStep  command,  except  that  the 
“Timestep”  is  the  “Tracking  time  per  timestep”  on  the  Advanced  Particle  Parameters  dialog 
box. 


Charge_Sur faces 
ReadOb j  ect 
OpenDatabase 

Se t Ini tialConduc tor Potentials 
Set Initial Potentials 
WritePotentials 
Embed_Ob j  ect_in_Grid 
Potent ials_in_Space 
Create_Particles 
Track_Particles 
Charge_Sur faces 

OpenDatabase 

InitializeCalculation 

SetEnvironment 

Environment 

type  LEO 

nel  densityvalue 
tel  temperaturevalue 
UseTrackedCur rents 
ZeroCurDerivAlgorithm 
PrepareChargeMatrix 
DoTrackTimeStep 
Potent ials_in_Space 
Track_Par  tides 
Charge_Sur faces 

OpenDatabase 

InitializeCalculation 

SetEnvironment 

Environment 

type  LEO 

nel  densityvalue 
tel  temperaturevalue 
UseTrackedCur rents 
ZeroCurDerivAlgorithm 
PrepareChargeMatrix 
DoTrackTimeStep 
Potent ials_in_Space 
Track  Particles 


The  plasma  density  and  temperature  appear  in  the  locations  indicated  by  italics.  The  value  of  the 
“UpdateTime”  keyword  in  the  “Track  Particles”  command  is  “No.” 
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The  tracked  current  density  appears  to  be  zero  on  the  3D  Results  tab  after  running  the  Charge 
Surfaces  module  because  the  tracked  current  is  associated  with  the  previous  timestep  and  only 
the  present  timestep  is  displayed.  The  charging  current  density  (which  for  this  calculation  is  the 
tracked  current  density)  is  associated  with  the  present  timestep  and  therefore  can  be  displayed. 
Contributions  to  the  net  current  from  photoemission  and  secondary  electrons  are  ignored.  Under 
some  conditions,  both  of  these  terms  can  be  important. 

A  number  of  minor  changes  were  made  to  the  way  the  histories  of  potentials  and  currents  are 
saved  in  order  to  avoid  the  display  of  possible  confusing  results.  The  two  most  visible  changes 
are  initialization  of  the  current  arrays  when  the  “SetlnitialPotentials”  command  is  executed  and 
the  saving  of  the  initial  surface  potentials  and  currents.  Results  for  time  “0.0”  are  now  displayed 
on  the  Results  tab. 

There  are  two  aspects  of  the  calculation  to  be  verified.  First,  that  the  current  to  surfaces 
computed  during  tracking  is  used  to  compute  the  change  in  potential  and,  second,  that  the  current 
is  being  used  correctly  in  the  calculation  of  the  change  in  surface  potential.  The  following 
calculations  were  used  to  verify  that  both  of  these  calculations  are  done  as  anticipated. 

4.1  Discharge  Conducting  Sphere  of  Known  Capacitance 

The  statement  of  the  test  case  is  as  follows: 

Object:  0.1-m  radius  “sphere”  shown  in  Figure  22. 

Environment:  1010m'3,  1  eV 
Initial  potential:  -100  V 

Space  potentials:  Hybrid  PIC  charge  density  model 

Particles:  Initialize  with  uniform  distribution  of  Hydrogen  ions,  5  x  10'7  s  timestep 
Charging:  Fifty  steps  of  above  script  without  Environment  specification. 
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Figure  22.  Sphere  Object  Used  in  Example 

The  current  collected  at  each  timestep  is  shown  in  Figure  23.  The  tracked  current  is  that  printed 
out  by  the  Track  Particles  module  and  the  charging  current  is  that  stored  by  the  Charge 
Surfaces  module.  This  figure  verifies  that  the  current  deposited  on  surfaces  by  the  Tracked 
Particles  module  is  used  by  the  Charge  Surfaces  module  to  compute  the  change  in  the  surface 
potentials. 


Figure  23.  Current  Collected  at  Each  Timestep  in  Discharge  Calculation 


The  change  in  potential  of  a  sphere  of  radius  r  in  vacuum  due  to  an  incident  current  I  during  a 

It 

timestep  of  length  r  is  given  by  Acjj  =  -7-  where  C  =  47ts0l'  .  Figure  24  compares  the 

potential  computed  by  Nascap-2k  with  the  potential  computed  from  the  capacitance  of  a  sphere 
of  radius  0.1  m  and  the  current  shown  in  Figure  23.  The  figure  also  shows  the  potential 
computed  for  the  capacitance  of  a  sphere  of  radius  0.09285  m. 
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Figure  24.  Potential  of  Sphere  During  Discharge 


The  ion  current  to  the  sphere  does  not  stop  immediately  when  the  sphere  goes  positive.  The  ions 
continue  to  move  toward  the  sphere,  and  to  be  collected,  due  to  the  particle  momentum.  When 
the  electric  field  is  high  enough  for  long  enough,  the  ions  are  all  moving  away  from  the  sphere 
and  the  current  goes  to  zero.  The  sphere  potential  stays  constant  after  this.  In  real  plasma,  the 
more  mobile  electrons  would  almost  completely  eliminate  the  overshoot  effect. 

The  positions  of  macroparticles  for  Z  values  between  0.0  m  and  0.03  m  at  5  ps,10  ps,  15  ps,  20 
ps,  and  25  ps,  are  shown  in  Figure  25  through  Figure  29.  In  the  first  figure,  the  positive  ions  can 
be  seen  moving  toward  negative  potential  sphere.  The  ions  continue  to  move  toward  the  sphere 
until  the  sphere  potential  is  positive  enough  long  enough  for  the  ions  to  begin  moving  away  from 
the  sphere. 


Figure  25.  Positions  of  Macroparticles  for  Z  Values  Between  0.0  and  0.03  m  at  5  ps. 
Particles  within  Sheath  Moving  Toward  Sphere 
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Figure  26.  Positions  of  Macroparticles  for  Z  Values  Between  0.0  and  0.03  m  at  10  ps 


Figure  27.  Positions  of  Macroparticles  for  Z  Values  Between  0.0  and  0.03  m  at  15  ps 
Surface  Potential  Near  Zero  -  Particles  Moving  Toward  Sphere 


Figure  28.  Positions  of  Macroparticles  for  Z  Values  Between  0.0  and  0.03  m  at  20  ps 
Surface  Potential  Positive  -  Current  Goes  to  Zero  as  Potentials  Overcome  Particle 

Momentum 
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Figure  29.  Positions  of  Macroparticles  for  Z  Values  Between  0.0  and  0.03  m  at  25  ps 
Particles  Moving  Away  from  Positive  Potential  Sphere 

The  results  of  repeating  the  above  calculation  with  the  inclusion  of  the  electron  current  are 
shown  in  Figure  30  and  Figure  31.  Without  the  contribution  of  the  electron  current,  the  sphere 
rises  to  nearly  10  V  positive  before  the  ion  current  drops  to  zero.  With  zero  incident  current,  the 
sphere  potential  remains  at  10  V  positive.  With  the  electron  current,  the  total  current  drops  to 
zero  as  the  potential  becomes  positive,  and  the  potential  is  held  near  zero  by  the  incident 
electrons.  The  ion  current  remains  fairly  constant  rather  than  dropping  to  zero,  as  the  ions  are  no 
longer  repelled.  Since  the  ion  current  is,  on  average,  slightly  less  than  the  raw  electron  thermal 
current,  the  potential,  on  average,  is  slightly  negative. 


Figure  30.  Comparison  of  Sphere  Potential  During  Discharge  with  and  without  an  Analytic 

Electron  Current  Included  in  the  Calculation 
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Figure  31.  Comparison  of  Current  Collected  at  Each  Timestep  in  Discharge  Calculation 
With  and  Without  an  Analytic  Electron  Current  Included  in  the  Calculation 

4.2  Differential  Charging  of  Insulating  Sphere  Surface 

The  statement  of  the  test  case  is  as  follows: 

Object:  0.1 -m  radius  Teflon  “sphere”  with  the  same  geometry  as  shown  in  Figure  22.  A 
dielectric  constant  of  2  and  a  thickness  of  0. 1  m  (completely  unphysical)  were  used  to  make 
the  capacitance  across  the  Teflon  near  the  capacitance  to  infinity. 

Environment:  1010  m'3,  1  eV 

Underlying  conductor  potential:  Fixed  at  -100  V 

Space  potentials:  Hybrid  PIC  charge  density  model 

Particles:  Initialize  with  uniform  distribution  of  Hydrogen  ions,  5  x  10'7  s  timestep 

Charging:  Thirty  steps  of  above  script  with  the  addition  of  a  “FixGroundPotential”  command 
with  an  argument  of  -100  V. 

The  current  collected  at  each  timestep  is  shown  in  Figure  32.  The  tracked  current  is  that  printed 
out  by  the  Track  Particles  module  and  the  charging  current  is  that  stored  by  the  Charge 
Surfaces  module.  A  surface  area  of  0.121  m2  is  used  to  make  these  currents  match. 
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Figure  32.  Current  Collected  at  Each  Timestep  in  Discharge  Calculation 

The  capacitance  between  the  surface  and  infinity  and  the  capacitance  between  the  surface  and  the 
underlying  conductor  are  in  parallel,  therefore,  the  capacitance  that  controls  the  charging  of  each 
Teflon  surface  is  their  sum.  The  capacitance  to  the  entire  surface  is  then  the  sum  of  the 
contributions  of  each  surface  and  equals  the  capacitance  between  the  sphere  and  infinity, 

C’  =  4nE0r  ,  and  the  capacitance  between  the  surface  and  the  underlying  conductor, 

C  =  KE0A/d ,  where  k  is  the  relative  dielectric  constant  of  the  Teflon,  A  is  the  surface  area, 
and  d  is  the  thickness  of  the  Teflon. 

Figure  33  compares  the  potential  computed  by  Nascap-2k  with  the  potential  computed  from  the 
ideal  capacitance  and  the  current  shown  in  Figure  32. 


Figure  33.  Potential  of  Sphere  During  Discharge 
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Figure  34.  Surface  Potentials  on  Insulating  Sphere  after  15  ps 
4.3  Further  Work 

Two  additional  features  to  the  analytic  electron  current  should  be  considered. 

1 .  Consider  the  electron  collection  for  positive  potentials.  Presently  it  is  enhanced  by  the 
three-dimensional  orbit  limited  (1+V/T).  However,  for  short  Debye  length  it  should  not 
be  enhanced.  Also,  there  are  cases  where  the  two-dimensional  formula  (something  like 
( 1+V/T) /z)  would  be  more  appropriate. 

2.  In  some  significant  cases  (e.g.,  VLF  antenna)  electron  collection  by  the  positive  surfaces 
is  blocked  by  the  negative  potential  sheath.  There  ought  to  be  a  way  to  approximate  this 
effect  using  the  electric  field  and  maybe  the  Debye  length,  but  it  is  not  obvious. 

The  inclusion  of  the  contributions  to  the  net  current  from  photoemission  and  secondary  electrons 
should  be  considered. 

5.  SELF-CONSISTENT  POTENTIALS  AND  CHARGE  DENSITIES  FOR  CHAWS 
PROBLEM  USING  HYBRID  PIC  APPROACH 

In  the  1990s  we  developed  a  technique  to  self-consistently  compute  the  potentials  and  charge 
densities  about  the  Wake  Shield  Facility,  a  free  flying  shuttle  payload,  with  the  biased  CHAWS 
probe.  The  technique,  referred  to  as  the  full  trajectory  approach,  is  to  iteratively  solve  for 
potentials  and  ion  charge  densities.  First  space  potentials  are  computed.  Then  a  thermal 
distribution  of  macroparticles  representing  ram  ions  is  tracked  from  the  problem  boundaries. 
Each  macroparticle  is  tracked  until  it  either  leaves  the  computational  space  or  hits  a  surface.  At 
each  step,  charge  is  deposited  into  the  volume  element  containing  the  macroparticle  at  the  end  of 
the  step.  After  all  of  the  macroparticles  are  tracked,  the  resulting  charge  density  is  used  to  solve 
Poisson’s  equation  for  the  space  potentials.  Ion  charge  densities  and  potentials  are  iteratively 
computed  until  a  steady  state  solution  is  reached.  To  reach  convergence,  we  found  that  at  each 
step  it  is  necessary  to  use  a  charge  density  composed  of  70%  of  the  previous  charge  density  and 
30%  of  the  most  recently  computed  charge  density.  A  similar  sharing  of  candidate  potential 
solutions  within  a  single  potential  computation  is  also  done.  The  potentials  and  charge  densities 
that  result  from  this  computation  are  shown  in  Figure  35  and  Figure  36. 
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Figure  35.  Potentials  in  Space  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  density  Computed  Using  Full  Trajectory  Approach 
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Figure  36.  Ion  Charge  Densities  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Full  Trajectory  Approach 
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Several  improvements  have  been  made  to  Nascap-2k,s  particle-in-cell  capabilities,  such  that  now 
the  same  problem  can  be  solved  using  a  Hybrid  PIC  approach.  The  original  implementation  of 
the  charge  stabilization  approach  for  the  solution  of  potentials  in  Hybrid  PIC  problems  was 
found  to  be  inadequate.  With  the  implementation  of  N2kDB,  we  are  able  to  save  both  nodal  and 
volume  charge  densities,  which  allows  us  to  compute  the  electron  charge  density,  including 
charge  stabilization,  for  a  Hybrid  PIC  problem  using  the  same  formulas  as  for  the  full  trajectory 
approach. 

Using  the  new  formulation  for  computing  the  electron  contribution  to  the  charge  density  with 
charge  stabilization,  we  found  that  while  the  solution  is  convergent,  the  tracked  ions  contribute 
little  to  the  charge  density  behind  the  WSF  disk.  This  is  because  the  method  originally  developed 
to  inject  macroparticles  at  the  problem  boundary  in  the  full  trajectory  approach  does  a  better  job 
capturing  the  wings  of  the  thermal  distribution  and  creates  a  higher  density  of  macroparticles 
near  the  disk  edge  than  does  the  method  originally  developed  for  Hybrid  PIC  problems. 

Another  improvement  to  the  standard  Hybrid  PIC  approach  is  “orbit  averaging.”  In  this 
approach,  rather  than  depositing  charge  to  the  nodes  for  the  ion  density  contribution  of  the 
charge  density  and  to  volume  elements  for  the  computation  of  the  electron  contribution  to  the 
charge  density  and  the  charge  stabilization  at  each  timestep,  charge  is  deposited  at  each  substep 
of  the  trajectory  and  accumulated  through  the  timestep.  This  allows  for  timesteps  longer  than  the 
time  to  cross  a  single  volume  element.  Using  orbit  averaging,  the  charge  density  behind  the  disk 
is  represented  a  little  better,  although  the  plasma  still  fills  the  wake  too  quickly  on  the  side  away 
from  the  probe.  In  addition,  the  potential  solution  is  less  well  behaved.  The  results  for  the  above 
approaches  were  reported  in  the  May- August  2009  quarterly  report,  with  plots  of  potential  and 
ion  charge  densities  illustrating  the  behavior  described  above. 

Recently,  the  flexible  macroparticle  specification  coding,  originally  developed  for  full  trajectory 
calculations,  was  modified  so  that  it  can  also  be  used  in  Hybrid  PIC  calculations.  Note  that 
presently,  the  new  thermal  distribution  splitting  capability  is  inconsistent  with  the  other  particle 
splitting  algorithms.  With  a  better  representation  of  the  thermal  distribution  and  a  spatial  particle 
distribution  with  increased  density  of  macroparticles  near  the  disk  edge,  the  potentials  and  charge 
densities  computed  using  the  Hybrid  PIC  approach  are  closer  to  those  obtained  using  the  full 
trajectory  approach.  The  results  are  shown  in  Figure  37  and  Figure  38,  with  noticeably  increased 
agreement  with  the  results  of  the  full-trajectory  approach  (see  Figure  35  and  Figure  36)  as 
compared  to  previous  efforts. 
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Figure  37.  Potentials  in  Space  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  5  ps 

Timestep,  after  1000  Timesteps  (5  ms) 
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Figure  38.  Ion  Charge  Densities  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  5  ps 

Timestep,  after  1000  Timesteps  (5  ms) 
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With  orbit  averaging,  the  results  of  this  static  problem  should  be  nearly  independent  of  the 
timestep  used.  With  a  spacecraft  velocity  of  7.8  km/sec,  it  takes  650  Q  for  an  ion  to  cross  the 
grid.  In  calculations  with  timesteps  near  and  larger  than  650  Q,  most  of  the  ions  cross  the  entire 
grid  within  each  timestep — increasing  the  resemblance  to  the  original  static  full  trajectory 
approach.  The  results  for  different  timesteps  are  shown  in  Figure  39  through  Figure  44.  As  the 
timestep  is  increased  the  wake  region  becomes  more  compact  and  the  results  more  unstable.  In 
the  shorter  timestep  calculations,  the  probe-side  potential  contour  shape  and  potential  gradient 
are  fairly  stable  from  timestep  to  timestep — after  the  initial  macroparticles  are  able  to  cross  the 
entire  gird.  With  longer  timesteps  the  results  are  less  stable.  The  potentials  are  less  smooth  and 
they  vary  more  from  timestep  to  timestep. 


Figure  39.  Potentials  in  Space  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  50  ps 

Timestep,  after  120  Timesteps  (6  ms) 
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Figure  40.  Ion  Charge  Densities  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  50  ps 

Timestep,  after  120  Timesteps  (6  ms) 
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Figure  41.  Potentials  in  Space  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  500 

ps  Timestep,  after  20  Timesteps  (10  ms) 
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Figure  42.  Ion  Charge  Densities  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  With  Orbit-averaging  Approach  with  a  500 

ps  Timestep,  after  20  Timesteps  (10  ms) 
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Figure  43.  Potentials  In  Space  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  5  ms 

Timestep,  after  10  Timesteps  (50  ms) 
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Figure  44.  Ion  Charge  Densities  for  CHAWS  Problem  -  Self-consistent  Potentials  and 
Charge  Density  Computed  Using  Hybrid  PIC  with  Orbit-averaging  Approach  with  a  5  ms 

Timestep,  after  10  Timesteps  (50  ms) 

Another  way  to  view  the  fidelity  of  the  simulation  is  the  current  to  the  CHAWS  probe.  The  probe 
current  from  the  full  trajectory  calculation  and  from  Hybrid  PIC  orbit-averaged  calculations  with 
5,  50,  200,  500,  1000,  and  5000  ps  timesteps  are  shown  in  Figure  45  as  a  function  of  timestep  (or 
iteration  for  full  trajectory).  With  a  timestep  of  650  ps  an  ion  crosses  the  computational  grid  in  a 
single  timestep,  and  so  these  calculations  span  static  and  dynamic  extremes  of  many  timesteps 
per  grid  crossing  to  several  grid  crossings  per  timestep. 

The  full  trajectory  approach  yields  a  rapidly  and  smoothly  converging  probe  current.  Note  that 
the  densities  are  computed  using  the  “sharing”  approach:  70%  of  the  previous  charge  density  and 
30%  of  the  most  recently  computed  charge  density  as  described  above.  A  few  calculations  were 
done  with  sharing  of  charge  densities  between  timesteps  in  the  same  manner  as  for  the  full 
trajectory  calculations.  The  probe  current  is  unaffected  by  the  presence/absence  of  sharing. 
(Those  with  sharing  denoted  ‘WS’  in  the  Figure  45  legend).  Figure  46  illustrates  the  convergence 
for  three  shorter- timestep  calculations  (5,  50  and  200  ps  timesteps),  plotting  probe  current  as  a 
function  of  simulation  time. 
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Figure  45.  Ion  Current  to  Probe  as  a  Function  of  Timestep,  for  Various  Timestep  Values 
and  as  a  Function  of  Iteration  for  the  Full  Trajectory  Calculation 


Figure  46.  Ion  Current  to  Probe  as  Function  of  Simulation  Time  for  5,  50,  and  200  ps 

Timesteps 
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As  seen  from  Table  1,  the  total  running  time  using  the  Hybrid  PIC  approach  increases 
significantly  when  a  large  number  of  timesteps  are  needed.  Note  that  large-timestep  calculations 
experience  greater  step-to-step  noise  even  at  late-time  convergence,  consistent  with  the  less 
stable  potentials  as  seen  in  Figure  39,  Figure  41,  and  Figure  43. 

Table  1.  Representative  Runtimes  for  Different  Timestep  Lengths 


Timestep/Params 

Minutes  to  Convergence 

Full  Traj 

160 

5  ps 

1800 

500  ps 

340 

5000  ps 

250 

6.  ORBIT  AVERAGING 

We  are  using  current  collection  by  a  sphere  to  explore  the  advantages  and  drawbacks  of  Nascap- 
2k’s  orbit  averaging  technique  to  compute  charge  densities  as  part  of  a  Hybrid  PIC  calculation  of 
potentials  in  space. 

We  computed  ion  collection  by  a  100  V,  0.1  m  radius  sphere  from  a  10  m'  ,  1  eV  plasma  in 
four  different  ways.  The  results  are  shown  in  Figure  47. 


Figure  47.  Ion  Current  to  -100  V  Sphere  Computed  Using  Four  Different  Techniques 
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The  equilibrium  current  is  computed  by  using  the  analytic  space  charge  approximation  and 
tracking  macroparticles  from  a  0.7  V  sheath  edge.  (Since  all  the  macroparticles  are  collected,  this 
result  is  just  the  area  of  the  sheath  surface.)  Since  this  approach  doesn’t  include  presheath 
enhancement,  the  result  is  then  multiplied  by  a  factor  of  1.45,  to  give  52.5  pA. 

The  Hybrid  PIC  solution  is  computed  by  initially  filling  the  computational  space  with  a  thermal 
distribution  of  macroparticles,  tracking  all  of  the  macroparticles  for  0.5  ps,  computing  potentials 
using  the  ion  charge  density  at  the  end  of  the  tracking  step,  and  then  repeating  the  tracking  and 
potential  computation  steps  for  0.7  ms  (1400  timesteps).  When  a  macroparticle  reaches  a  grid 
boundary,  it  is  split  into  multiple  macroparticles,  if  it  carries  too  much  charge  for  the  new  grid.  A 
thermal  distribution  of  new  macroparticles  are  injected  from  the  boundary  every  tenth  timestep 
in  order  to  compensate  for  the  current  collected  by  the  sphere.  The  equilibrium  solution  is 
reached  in  about  0.2  ms.  This  is  about  as  long  as  needed  for  most  of  the  macroparticles  incident 
on  the  boundary  at  the  first  timestep  to  reach  the  sphere. 

The  Hybrid  PIC,  orbit-averaged  solution  is  computed  in  a  similar  manner  to  the  Hybrid  PIC 
solution.  The  difference  is  that  rather  than  depositing  charge  to  the  nodes  for  the  ion  density 
contribution  of  the  charge  density  and  to  volume  elements  for  the  computation  of  the  electron 
contribution  to  the  charge  density  and  the  charge  stabilization  at  each  PIC  timestep,  charge  is 
deposited  at  each  step  of  the  trajectory  and  accumulated  through  the  timestep.  This  allows  for 
timesteps  longer  than  the  time  to  cross  a  single  volume  element.  With  a  5  ps  timestep  and 
boundary  current  injection  at  every  timestep,  the  current  to  the  sphere  is  similar  to  that  given  by 
the  standard  Hybrid  PIC  approach.  As  macroparticles  are  injected  every  timestep,  rather  than 
every  10  timesteps,  the  number  of  macroparticles  is  similar  to  that  in  the  standard  Hybrid  PIC. 
The  total  running  time  is  about  the  same,  as  the  potential  computation  is  only  5%  to  10%  of  the 
total  execution  time. 

We  also  computed  the  current  to  the  sphere  using  the  Hybrid  PIC,  orbit-averaged  approach  and  a 
50  ps  timestep.  This  calculation  has  large  variations  in  the  collected  current  from  cycle  to  cycle. 
Using  the  present  algorithms,  the  Hybrid  PIC  solution  to  ion  current  collection  by  a  negative 
sphere  at  constant  potential  reaches  the  equilibrium  solution  rapidly.  The  orbit  averaging 
technique  gives  the  same  results  and  requires  about  the  same  amount  of  computational  time  if  the 
integration  time  is  a  fraction  of  the  time  need  for  the  bulk  of  the  ions  to  either  be  collected  or 
leave  the  computational  space. 

The  role  of  initialization  and  adjustments  to  the  technique  needed  to  accommodate  timesteps 
comparable  to  and  longer  than  the  time  needed  to  cross  the  computational  space  remain  to  be 
explored. 


7.  PSEUDOPOTENTIAL  ALGORITHMS  FOR  SIMULATION  OF  VLF  PLASMA 
ANTENNA  CURRENT  FLOW 

We  developed  pseudopotential  approaches  (similar  to  the  velocity  potential  used  to  describe 
potential  flow  [1]  in  fluid  dynamics)  to  the  computation  of  transverse  surface  currents  on  a 
spacecraft  acting  as  an  antenna  and  to  the  computation  of  electron  currents  in  the  sheath  and  near 
sheath  about  the  spacecraft.  These  approaches  have  been  implemented  within  Nascap-2k. 
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The  transverse  surface  current  algorithm  in  Nascap-2k  is  used  to  calculate  the  current  flowing 
along  an  antenna  element  or  other  object,  accounting  correctly  for  both  local  capacitance  and 
incident  plasma  current.  The  surface  is  divided  into  many  surface  elements.  The  change  in 
charge  on  a  surface  element  during  a  timestep  is  that  which  is  needed  to  accomplish  the  change 
in  surface  electric  field,  which  is  obtained  from  the  Nascap-2k  potential  calculation.  The  current 
to  each  surface  element  consists  of  the  transverse  surface  current  from  neighboring  surface 
elements  and  the  current  provided  by  the  plasma,  which  is  provided  by  the  Nascap-2k  PIC 
(particle-in-cell)  simulation.  The  current  continuity  equation  is  solved  using  a  pseudopotential 
approach.  As  a  boundary  condition,  one  surface  element  of  each  antenna  element  is  specified  as 
connected  to  the  biasing  power  supply.  With  pseudopotential  values  assigned  to  each  surface 
element,  the  solution  provides  the  currents  crossing  edges  between  surfaces  needed  to  satisfy  the 
problem  dynamics.  This  approach  excludes  solutions  with  circulating  currents.  The  vector 
transverse  surface  current  in  each  surface  element  is  taken  as  that  which  provides  the  best  fit  to 
the  edge  currents. 

The  volume  current  algorithm  is  used  to  calculate  electron  currents  within  and  near  the  sheath 
about  a  VLF  antenna  or  other  high-voltage  object.  Within  the  Nascap-2k  framework,  local 
equilibrium  electron  densities  are  generated  for  each  volume  element  as  part  of  the  PIC  ion 
dynamics  algorithm,  and  their  time  derivatives  are  the  main  drivers  for  volume  currents.  Space 
outside  the  calculation  boundary  can  act  as  either  a  source  or  a  sink  for  electrons,  and  object 
surfaces  may  act  as  electron  sinks. 

The  example  used  here  is  a  Nascap-2k  simulation  of  the  DSX  transmitting  antenna  in  low 
density  hydrogen  plasma.  Parameters  of  the  calculation  are  shown  in 

Table  2.  Note  that  the  1000  V  square  wave  bias  is  approximated  by  the  sum  of  the  first  and  third 
harmonics,  so  that  the  actual  peak-to-peak  voltage  is  closer  to  1200  V. 

Table  2.  Parameters  For  DSX  Antenna  Simulation 


Antenna  Length  (tip  to  tip) 

81 

Antenna  Diameter 

0.1  m 

Bias  Voltage  (Peak  to  Peak) 

1000  V 

Bias  Frequency 

10  kHz 

Plasma  Density 

o 

oo 

3, 

Plasma  Temperature 

1  eV 

7.1  Mathematical  Formulation 


The  basic  equation  to  be  solved  is  V  •  J  +  ^  =  0,  where  J  is  the  surface  or  volume  current  density 

due  to  electrons,  and  p  is  the  surface  or  volume  electron  charge  density.  Note  that  the  solution  for 
the  current  is  non-unique  to  the  addition  of  any  divergence-free  current  field.  We  assume  that, 
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provided  appropriate  boundary  conditions  are  implemented,  such  circulating  currents  are  not  of 
concern  to  the  problem  being  solved. 

A  solution  having  no  circulating  currents  is  guaranteed  if  we  assign  a  pseudopotential  value  to 
each  element  of  surface  or  volume  and  take  the  current  across  their  common  boundary  (edge  or 
face)  to  be  proportional  to  the  difference  in  their  pseudopotentials.  (In  that  case,  current  that 
flows  “downhill”  would  need  to  flow  back  “uphill”  in  order  to  “circulate.”)  The  current  is  also 
taken  as  proportional  to  the  interface  area  (edge  length  or  face  area)  and  inversely  proportional  to 
the  distance  between  the  centers  of  the  elements. 

The  result  of  the  above  treatment  is  a  matrix  that  multiplies  the  vector  of  pseudopotentials  to 
describe  the  buildup  of  charge  in  the  surface  or  volume  elements,  which  is  then  set  equal  to  the 
source  vector  generated  from  Nascap-2k  results  at  each  timestep.  This  system  of  equations  is 
then  solved  using  the  ICCG  (Incomplete  Cholesky  Conjugate  Gradient)  [2]  algorithm. 

7.2  Application  to  Surface  Currents 

The  spacecraft  has  two  or  more  antenna  elements  (represented  by  conductor  numbers  in 
Nascap-2k )  that  are  biased  in  a  programmed  way  (generally  periodic)  by  an  amplifier.  Typically 
the  spacecraft  will  have  passive  elements  as  well.  The  amplifier  injects  current  at  a  specified 
location  on  each  antenna  element  to  maintain  the  programmed  voltage.  The  injected  current, 
modified  by  current  collected  from  the  plasma,  flows  on  the  spacecraft  surfaces  so  as  to  maintain 
the  specified  potential  difference  between  each  surface  element  and  ambient  plasma. 

Figure  48  shows  graphically  the  quantities  that  make  up  the  surface  current  equation  for  each 
surface  element.  Moving  the  surface  currents  to  the  left  side  of  the  equation,  and  all  other 
(known)  quantities  to  the  right,  and  expressing  the  edge  currents  as  proportional  to 
pseudopotential  differences,  results  in  a  sparse  symmetric  matrix  that  is  amenable  to  solution 
using  the  ICCG  algorithm  once  the  surface  element  connected  to  the  power  supply  (denoted  the 
“injection  element”)  is  set  to  a  fixed  potential.  The  current  required  to  bias  the  conductor  may  be 
obtained  by  evaluating  the  equation  associated  with  the  injection  element  and  can  be  verified  to 
be  the  sum  of  the  total  change  in  charge  less  plasma  currents  for  all  the  remaining  surfaces  of  the 
conductor.  Surface  currents  are  solved  separately  for  each  electrically  isolated  component. 


Total  change  in  charge  Plasma  Currents  Surface  Currents 


Figure  48.  The  Change  in  Charge  on  a  Surface  Element  is  Made  Up  of  Plasma  Currents 

and  Surface  Currents 
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The  following  figures  show  surface  currents  from  a  simulation  of  the  DSX  model  at  the  time  of 
peak  surface  current  in  the  antenna  elements,  i.e.,  the  surface  current  while  the  bias  polarity  is 
switching.  A  close-up  of  the  portion  of  the  antenna  elements  near  the  roots  with  surfaces  color- 
coded  by  surface  current  magnitude  is  shown  in  Figure  49.  The  height  of  each  cone  is 
proportional  to  the  square  root  of  the  surface  current  density  and  to  the  square  root  of  the  surface 
element  area.  Note  that  the  current  flows  in  the  same  direction  on  both  antenna  elements.  Figure 
50  shows  the  full  length  of  both  antenna  elements.  The  current  is  a  maximum  at  the  antenna 
element  roots  (i.e.,  at  the  power  supply)  and  decreases  toward  the  tips.  The  current  magnitude  as 
a  function  of  location  along  the  antenna  element  is  plotted  in  Figure  51,  with  the  spatial 
derivative  of  surface  current  plotted  using  the  scale  on  the  right.  This  shows  that  the  current 
varies  linearly  over  most  of  the  antenna  elements,  with  a  large  amount  of  capacitive  loading  near 
the  root  and  a  smaller  amount  near  the  antenna  element  tip.  Figure  52  illustrates  the  time 
dependence  of  the  surface  current,  plotted  together  with  the  applied  voltage.  As  expected  for  a 
capacitive  system,  the  surface  current  is  the  time  derivative  of  the  potential  (with  very  small 
modification  due  to  incident  plasma  current),  and  this  is  largest  when  the  two  antenna  elements 
are  switching  polarity,  and  strongly  emphasizes  the  third  harmonic  frequency. 
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Figure  49.  Graphical  Display  of  Transverse  Surface  Currents,  with  Equation  Relating 
Element  Surface  Currents  to  Vector  Potential 
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Figure  50.  Another  Graphical  Display  of  Surface  Current  Flow  on  Antenna  Elements, 
Showing  Current  Maximum  at  Antenna  Element  Roots  Decreasing  Toward  Tips 
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Figure  51.  Plot  of  Surface  Currents  Corresponding  to  One  Antenna  Element  of  Figure  50, 
Showing  that  Current  Varies  Nearly  Linearly  Along  the  Antenna  Element,  With 

Capacitive  Loading  at  Ends 
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Figure  52.  Variation  of  Average  Transverse  Current  in  One  Antenna  Element  Over  Two 

Full  Cycles 
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7.3  Application  to  Volume  Electron  Currents 

Generation  of  propagating  waves  in  plasma  is  due  primarily  to  electron  currents  flowing  in  the 
bulk  cold  plasma.  While  plasma  ion  currents  can  be  generated  directly  by  the  PIC  ion  simulation 
used  to  obtain  the  sheath  dynamics,  it  is  impractical  to  simulate  electron  macroparticles  to  obtain 
electron  currents.  The  calculation  we  use  as  an  example  is  a  simulation  using  greater  than  6><1 07 
ion  macroparticles.  If  electrons  were  tracked  as  well  (instead  of  using  a  “fluid”  or  “local 
equilibrium”  approximation),  a  much  larger  number  of  electron  macroparticles  would  be 
required,  as  well  as  a  shorter  timestep.  The  plasma  ion  current  patterns,  obtained  directly  from 
the  ion  tracking,  are  discussed  in  Section  10. 

The  main  requirement  for  application  of  the  pseudopotential  approach  to  the  calculation  of 
volume  electron  currents  is  that  the  electron  density  for  each  volume  element  (as  used  in  the 
calculation  of  electrostatic  potential  during  the  PIC  calculation)  be  stored  (or  computable)  for 
each  timestep.  In  addition,  we  need  to  establish  boundary  conditions  at  the  external  boundary  and 
near  the  spacecraft,  develop  a  local  conductivity  tensor  that  reflects  the  non-uniform  and 
anisotropic  nature  of  the  plasma,  and  do  extra  bookkeeping  for  currents  crossing  grid  boundaries. 
We  have  thus  far  restricted  application  of  the  pseudopotential  method  to  cubic  (empty)  volume 
elements  and  allow  current  flow  through  the  six  faces  of  each  cube,  i.e.,  in  the  X,  Y,  and  Z 
coordinate  directions. 

The  volume  external  to  the  computational  space  is  assigned  a  single  pseudopotential  with  value 
fixed  at  zero.  This  allows  the  external  space  to  be  a  source  or  sink  of  electrons  at  each  element  of 
the  external  boundary  surface. 

In  Nascap-2k ,  the  object  is  surrounded  by  “special,”  (i.e.,  non-cubic)  volume  elements.  At  the 
interface  between  an  empty  element  and  a  special  element,  we  allow  the  special  element  to  be  a 
sink,  but  not  a  source,  of  electrons.  Specifically,  if  the  electric  field  in  the  empty  element  points 
toward  the  special  element,  we  take  the  electron  current  crossing  the  interface  to  be  zero.  If  the 
electric  field  points  away  from  the  special  element,  we  take  the  electron  current  crossing  the 
interface  to  be  the  plasma  thermal  current  at  the  empty  element’s  density.  This  treatment 
represents  the  collection  of  electrons  by  positive  surfaces. 

The  treatment  must  also  recognize  that  electron  current  flows  along  paths  with  an  electron 
population  that  can  support  such  current  and  not  along  paths  that  are  devoid  of  electrons  (such  as 
the  interior  of  a  sheath).  To  that  end,  we  set  the  conductivity  to  be  proportional  to  the  square  root 
of  local  electron  density  (i.e.,  to  the  electron  density  divided  by  a  scattering  rate  taken  to  be 
proportional  to  the  electron  plasma  frequency).  The  conductivity  is  expressed  as  a  tensor  so  that 
anisotropy  due  to  magnetic  field  can  be  taken  into  account. 

To  illustrate  the  behavior  of  this  algorithm,  we  use  a  DSX  transmitting  antenna  simulation 
having  the  parameters  shown  in  Table  2.  The  current  pattern  at  the  time  that  the  antenna 
elements  are  switching  is  shown  in  Figure  53.  At  this  point  in  time  both  antenna  elements  are  at 
-577  V,  with  the  upper  element  becoming  more  negative  so  that  its  sheath  is  growing  and 
expelling  electrons,  and  the  lower  positive  element  becoming  less  negative  so  that  its  sheath  is 
shrinking  with  electrons  filling  in  the  formerly  negative  space.  In  the  left  (X-component)  portion 
of  the  figure,  the  light  area  at  upper  left  indicates  current  moving  to  the  right,  or  electrons  to  the 
left,  away  from  the  upper  antenna  element.  Similarly,  the  light  and  dark  areas  in  the  remaining 
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quadrants  all  show  electrons  moving  away  from  the  upper  antenna  element  and  toward  the  lower 
antenna  element.  The  Y-component  currents  at  the  same  time  show  upward  current  (downward 
electron  motion)  in  the  neighborhood  of  the  sheath  edge.  Note  that  there  is  negligible  current  in 
the  interior  of  the  sheath.  If  the  conductivity  tensor  does  not  account  for  the  low  electron  density 
in  the  sheath,  the  entire  sheath  is  filled  with  substantial  upward  current,  which  is  unphysical 
since  electrons  are  excluded  from  the  region. 
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Figure  53.  X  (Left)  and  Y  (Right)  Component  of  Current  (Am'2)  when  Antenna  Elements 

are  Switching  -  See  Discussion  in  Text 

Figure  54  shows  the  electron  currents  generated  when  the  lower  antenna  element  is  positive  and 
collecting  electron  current.  The  outward  current  (inward  electron  flow),  which  is  needed  both  to 
supply  the  electrons  collected  by  the  antenna  element  and  to  fill  in  the  contracting  sheath,  fills  a 
large  volume  around  the  antenna.  When  the  currents  are  plotted  on  a  log  scale,  the  smaller 
outward  electron  flow  can  be  seen  near  the  still-expanding  sheath  edge  around  the  negative 
antenna  element.  Figure  55  shows  that  a  short  time  later  the  electron  current  has  reversed  sign  as 
the  lower  antenna  element  has  become  negative  and  is  expelling  electrons  from  its  sheath. 
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Figure  54.  X  Component  of  Electron  Current  on  Linear  and  Log  Scales  when  Lower 
Antenna  Element  is  at  +0.16  V  -  Electrons  are  Collected  by  Lower  Antenna  Element 
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Figure  55.  X  Component  of  Electron  Current  when  Lower  Element  Switches  from  Positive 
to  Negative  Potential  (Here  at  -13  V)  -  Electrons  are  Expelled  from  the  Near-Antenna 
Region  to  Form  a  Sheath,  and  Flow  Outward 
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Figure  53  (left  panel)  indicates  electrons  moving  outward  into  ambient  plasma  in  the  upper  half 
of  the  figure  and  inward  in  the  lower  half.  These  counterstreaming  electron  currents  are 
commonly  associated  with  whistler  waves.  The  right  panel  of  Figure  53  suggests  that  the  closing 
Y-component  current  is  small  at  the  boundary  so  that  the  counterstreaming  currents  have 
considerable  extent.  If  the  magnetic  field  is  in  the  X-direction,  then  the  closing  Y-component 
current  is  suppressed  and  the  extent  of  the  counterstreaming  currents  enhanced. 

7.4  Summary 

We  have  demonstrated  the  use  of  pseudopotential  methods  to  calculate  surface  currents  and 
volume  electron  currents.  The  source  terms  for  these  currents  are  obtained  from  Nascap-2k 
dynamic  sheath  calculations.  Boundary  conditions  are  imposed  based  on  physical  considerations. 

The  source  terms  for  the  surface  current  algorithm  are  the  rate  of  change  of  surface  electric  field 
at  each  surface  element  and  the  incident  plasma  currents.  The  fixed  pseudopotential  boundary 
condition  of  an  active  antenna  element  is  applied  at  the  location  of  the  power  supply  connection 
(antenna  element  root).  For  a  passive  object  component  any  surface  element’s  pseudopotential 
may  be  fixed,  as  this  object  component’s  total  charge  changes  only  due  to  incident  plasma 
currents. 

The  source  term  for  the  volume  electron  current  algorithm  is  the  calculated  rate  of  change  of 
electron  charge  in  each  volume  element.  The  external  space  is  held  at  fixed  pseudopotential  so 
that  it  can  act  as  either  a  source  or  a  sink  for  electrons.  Positive  surfaces  act  as  sinks  for 
electrons.  A  local  conductivity  tensor  is  used  to  reflect  the  ability  of  the  electron  population  to 
carry  current,  embodying  the  effects  of  magnetic  field  and  local  electron  density. 

A  surface  current  example  was  presented  showing  calculation  of  current  on  the  DSX  transmitting 
antenna.  Results  showed  that  the  current  varied  linearly  over  most  of  the  antenna  element,  with 
capacitive  loading  at  the  root  and  tip. 

Volume  electron  currents  were  illustrated  by  the  same  DSX  calculation  and  showed 
counterstreaming  electron  currents  in  the  distant  plasma. 


8.  PROTOTYPE  OF  NASCAP-2K  REALTIME 

Nascap-2k  RealTime  Prototype  computes  surface  potentials  on  spacecraft  in  response  to  tabular 
spectra  provided  in  real  time.  It  is  a  stand-alone  Java  application  and  uses  a  robust  version  of  the 
Boundary  Element  Method  (BEM)  charging  algorithms  developed  for  Nascap-2k  and  originally 
implemented  in  Java  in  the  SEE  Interactive  Spacecraft  Charging  Handbook.  The  charging 
algorithms  used  are  only  appropriate  to  the  low  density  plasma  environments  found  at 
geostationary  altitude.  The  algorithm  used  to  determine  the  sun  direction  is  also  only  appropriate 
to  spacecraft  at  geostationary  altitude. 

N as cap-2k  RealTime  Prototype  executes  in  response  to  the  command  line  prompt 

java  -jar  Nascap2kRT.jar  -prefix  MiniDSCS  -envFileName  SW199610500900sat_01.0.MSM  - 

satConfFileName  sateljCOl.Ol  where  the  options  are  specified  in  Table  3. 
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Nascap-2k  RealTime  Prototype  requires  four  input  files: 

.  Object  Toolkit  description  of  spacecraft,  including  specification  of  rotating  solar  arrays. 
{prefix  Obj  ect.xml) 

.  List  of  times,  plasma  environments,  and  positions  of  the  spacecraft.  The  environment  can 
be  expressed  either  as  a  table  of  fluxes  in  various  energy  bins  or  as  a  Maxwellian.  (MSM 
fde  or  XML  file) 

.  Previously  computed  potentials  for  each  surface  at  each  timestep  and  the  environments 
used  to  compute  the  potentials.  (pre/zxSteps.xml) 

•  Satellite  configuration  file  used  only  to  determine  the  heading  of  the  output  file. 

Nascap-2k  RealTime  Prototype  creates  two  output  files: 

•  Computed  potentials  for  each  surface  at  each  timestep  and  environment  used  to  compute 
the  potentials.  (pre/zxSteps.xml) 

•  Table  of  chassis  potential,  minimum  differential  potential,  and  maximum  differential 
potential  for  the  times  specified  in  the  input  file.  (SUR  file) 

The  format  of  the  MSM,  SUR,  and  configuration  files  are  given  in  Interface  Control  Document 
for  The  SEEFS  Satellite  Charging/Discharging  Product  dated  May  11,  2005.  The  format  of  the 
Maxwellian  environment  and  /ire/ixSteps.xml  files  are  given  in  the  examples. 

The  calculation  consists  of  five  steps: 

.  Object  initialization 

•  Read  environments 

.  Initialize  calculation 

•  Timestepping 

•  Write  output  files 


57 


Table  3.  Options  for  Running  Nascap-2k  RealTime  Prototype 


Option 

Meaning 

Default 

-prefix  prefix 

The  Object  ToolKit  file  for  the  object  must  be  in  the  run 
directory  as  prefix  Object.xml. 

Calculated  potentials  are  written  to  and  read  from 
prejixS  tep  s .  xml . 

None 

-envFileName  EnvFile 

Name  of  the  environment  file  to  use.  Either  an  MSM 
file  or  a  Maxwellian  in  an  xml  format. 

None 

-satConfFileName 

satConfFile 

Name  of  the  satellite  configuration  file  to  use  (e.g. 
satelCOl.Ol). 

None 

-dir  Directory 

Specifies  the  directory  in  which  all  input  and  output 
files  appear. 

Local 

directory 

-maxDataSaveTime 

maxSaveTime 

The  maximum  time  span  to  save  data  in  pre/ixSteps.xml 
file;  specified  in  hours. 

6  hours 

-timeStep  TStep 

Timestep  (in  seconds).  Presently  must  be  an  integer. 

10  seconds 

-advanceTime  advTime 

Minutes  to  continue  the  calculation  past  the  last 
environment  in  the  environment  file.  Presently  must  be 
an  integer. 

0  minutes 

-minsToUpConf 

confUpdateTime 

Number  of  minutes  the  calculation  runs  without  a  bad 
result  before  the  confidence  level  is  increased.  This  is  a 
real  number. 

1  minute 

8.1  Calculational  Steps 

8.1.1  Read  Command  Line  Input  and  Initialize  Object 

The  first  step  of  the  calculation  is  to  read  the  command  line  options  and  attempt  to  create  a  lock 
file,  prefix. Ick.  If  the  file  prefix.lck  already  exists,  the  code  exits.  This  test  insures  that  for  a 
single  directory  and  single  prefix,  only  one  instance  of  Nascap-2k  RealTime  Prototype  can 
execute  at  one  time.  The  code  then  reads  the  configuration  file  and  then  the  geometry  from  the 
Object  ToolKit  file  prefix Object.xml.  If  the  geometry  has  rotating  parts  like  solar  arrays,  the 
rotating  surfaces  are  specified  in  the  object  file  in  the  standard  way.  The  BEM  matrix  elements 
are  computed  for  the  object  without  rotation.  If  the  object  definition  includes  rotating 
components,  the  matrices  are  later  recomputed  for  each  sun  direction  provided. 

8.1.2  Read  Environments 

The  environments  are  read  from  the  file  specified  in  the  command  line.  That  the  environment  is 
given  as  a  table  of  fluxes  in  a  set  of  energy  bins  is  indicated  by  any  extension  except  xml  in  the 
file  name.  That  the  environment  is  specified  as  a  Maxwellian  is  indicated  by  an  extension  of  xml 
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in  the  file  name.  The  environment  is  then  specified  by  a  list  of  densities,  temperatures,  spacecraft 
positions,  and  associated  times. 

If  there  are  lines  in  an  MSM  file  that  do  not  produce  valid  environments,  these  lines  are  ignored 
and  the  code  uses  only  the  valid  lines.  Invalid  environments  can  occur  if  the  data  is  incomplete, 
has  negative  fluxes,  etc. 

The  spacecraft  position  is  used  to  determine  both  the  orientation  of  any  rotating  components  and 
the  direction  to  the  sun  for  the  computation  of  photoemission. 

Once  all  the  environments  are  read  from  the  file,  the  environments  are  ordered  by  time  for  use  in 
the  calculation. 

8.1.3  Initialize  Calculation 

The  initialization  step  begins  by  reading  any  previously  computed  results  from  the 
y>re//xSteps.xml  file. 

If  there  are  no  previous  results,  the  time  is  set  to  the  earliest  time  in  the  environment  file,  the 
environment  is  set  to  the  environment  at  that  time,  and  all  the  surface  potentials  are  set  to  0.0  V. 

If  there  are  previous  results,  the  environments  in  the  file  are  compared  with  those  used  in  the 
previous  calculations  for  the  corresponding  time.  For  tabular  environments,  the  comparison  is 
between  the  corresponding  flux  values,  bin  edge  values,  and  eclipse  indicators.  For  Maxwellian 
environments,  the  comparison  is  between  the  corresponding  densities  and  the  temperatures. 
Previously  computed  results  for  which  the  new  environment  is  different  from  the  one  used  in  the 
previous  computation  are  discarded.  The  time  is  set  to  the  latest  time  for  which  the  previously 
computed  results  are  kept.  The  potentials  on  all  the  surfaces  are  set  to  the  previously  computed 
values  at  this  time.  The  environment  is  set  to  the  new  environment  for  this  time  if  one  is 
available.  If  not,  the  environment  is  set  to  the  environment  for  this  time  saved  with  the  previous 
results. 

Finally,  any  rotating  surfaces  are  rotated  and  new  BEM  matrix  elements  are  computed. 

8.1.4  Time  Stepping 

The  heart  of  the  calculation  is  the  time  stepping  through  the  time  period  specified  in  the 
environment  file.  The  length  of  each  timestep  is  set  to  the  shorter  of  the  timestep  specified  in  the 
argument  list  and  the  time  between  the  present  time  and  the  time  associated  with  the  next 
specified  environment. 

At  each  timestep,  the  code  records  the  time,  the  maximum  potential,  the  minimum  potential,  the 
conductor  potentials,  the  confidence  level,  and  the  potentials  on  each  surface  in  xml  format.  This 
information  is  written  to  the  prefixSteps.xml  file  when  the  code  is  done  (not  during  execution). 
The  code  continues  to  run  until  the  present  time  reaches  the  time  associated  with  the  latest 
environment  in  the  input  file  plus  the  time  period  specified  in  the  argument  list. 
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At  the  end  of  each  timestep,  the  code  checks  for  a  new  environment  for  the  present  time.  If  one 
exists,  the  spacecraft  position  is  used  to  compute  the  position  of  any  rotating  surfaces  and  the  sun 
direction.  New  BEM  matrix  elements  are  computed  for  the  new  spacecraft  configuration. 

8. 1.4.1  Algorithm  for  the  Computation  of  Confidence  Level 

The  code  uses  four  confidence  levels;  0,  30,  50,  and  70.  When  a  new  calculation  is  started  and 
the  potentials  are  all  set  to  0,  the  confidence  level  is  set  to  0. 

If  the  potential  on  any  surface  after  a  timestep  is  NaN  (not  a  number),  all  potentials  are  set  to  0 
and  the  confidence  level  is  reset  to  0.  (We  expect  this  to  be  infrequent.) 

At  the  end  of  each  timestep,  if  there  is  a  surface  potential  >  20  V  or  <  -10,000  V,  the  confidence 
is  decreased  one  level.  If  the  chassis  potential  changes  more  than  1  kV  over  one  minute  while  the 
satellite  is  sunlit,  the  confidence  level  is  decreased  by  one.  If  the  differential  potential  increases 
or  decreases  by  more  than  500  V  over  one  minute,  the  confidence  level  is  decreased  by  one. 

If  the  confidence  level  has  not  changed  for  the  specified  period  of  time,  it  increases  by  one  level. 

8.1.4.2  Surface  Currents 


For  environments  specified  as  a  table,  the  surface  currents  are  computed  from  the  following 
expression. 
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The  first  term  corresponds  to  the  orbit  limited  collection  and  the  second  to  the  orbit  limited 
collection  from  energy  bins  that  are  partially  excluded. 


8.1.5  Create  Output  Files 

At  the  completion  of  execution,  the  code  writes  two  files;  the  pre/ixSteps.xml,  which  contains 
the  time,  the  maximum  potential,  the  minimum  potential,  the  conductor  potentials,  the 
confidence  level,  and  the  potentials  on  each  surface  at  each  timestep  in  XML  format,  and  the 
.SUR  file,  which  contains  the  chassis  and  differential  potentials  along  with  the  confidence 
estimate  at  the  times  specified  in  the  input  MSM  file. 


60 


8.2  Objects 


Three  geometric  objects  are  included;  a  Kapton  coated  sphere  (Figure  56),  a  DSCS-like 
spacecraft  with  only  insulating  surfaces  (Figure  57),  and  a  DSCS-like  spacecraft  with  ITO  coated 
solar  arrays  (Figure  58). 


Figure  56.  Spherically  Shaped  Object  Available  for  Calculations 


Figure  57.  DSCS-Like  Spacecraft  Geometry  Available  for  Calculations  -  the  Solar  Arrays 
Rotate  About  the  Long  Axis  in  Order  to  Track  the  Sun 
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Figure  58.  Second  DSCS-Like  Spacecraft  Geometry  Available  for  Calculations  -  the  Solar 
Arrays  Rotate  About  the  Long  Axis  in  Order  to  Track  the  Sun 

8.3  Supporting  Software 

For  convenience  in  viewing  results  or  plotting,  there  is  a  MakeDataFile.jar  that  takes  an  output 
xml  file  as  input  and  produces  a  text  file  in  which  each  line  has  the  time  in  seconds  from  the  first 
time,  the  chassis  potential,  the  minimum  potential,  and  the  maximum  potential. 

java  -jar  MakeDataFile.jar  prefix Steps.xml 

8.4  Verification 

In  order  to  verify  that  Nascap-2k  RealTime  Prototype  is  calculating  spacecraft  surface  charging 
correctly,  we  compared  results  computed  using  the  new  code  with  results  computed  using  the 
SEE  Interactive  Spacecraft  Charging  Handbook.  As  the  SEE  Handbook  has  restrictive  geometry, 
the  Object  ToolKit  object  developed  for  early  cross  code  comparisons  was  used.  The  spacecraft 
is  taken  to  be  at  0°  longitude  and  0°  latitude.  The  date  is  January  1,  2000.  The  environment  is  the 
NASA  Worst  Case  surface  charging  environment.  The  results  for  midnight  and  for  6  a.m.  (solar 
array  orientation  and  sun  direction)  computed  using  the  SEE  Handbook,  using  Nascap-2k 
RealTime  Prototype  and  a  Maxwellian  environment  expressed  as  a  density  and  temperature,  and 
using  Nascap-2k  RealTime  Prototype  and  a  Maxwellian  environment  expressed  as  a  table  of 
fluxes  and  energy  bins  are  shown  in  Figure  59  and  Figure  60.  The  results  can  only  be  as  close  as 
shown  if  all  aspects  of  the  charging  calculation  are  done  correctly.  The  same  comparison,  redone 
with  ten-second  timesteps  in  Nascap-2k  RealTime  Prototype,  is  shown  in  Figure  61  and  Figure 
62. 
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Potential  I  Vo  Its 


Midnight:  January  1,  2000^  NASA  Worst  Case  environment;  Using  Handbook  Timesteps 


■  Handbook-Chassis 
— * —  Hand  book -Min 
Handbook-Ma* 
Maxwe-lian-Chass  s 
— Mswe  lian-Min 
— ■ —  Ma*walian-Max 
—l — Table-Chassis 
^^Table-Min 
Table-Max 


Figure  59.  Comparison  of  Minimum,  Maximum,  and  Chassis  Potentials  as  a  Function  of 
Time  Computed  in  Three  Different  Ways  for  Midnight  on  January  1, 2000 


6  AM  :  January  1,2000:  NASA  Worst  Case  environment;  Using  Handbook  Timesteps 


-  Handbook-Chass 'S 

-  Handbook-Min 
Handbook-Ma* 
Mascwe'  ianChass  e 

-  Maawe&ian-Mi?’- 
-MaxweSian-Max 
-TaUfeOHSEis 
-Tafcte-Min 

Tto-Ma* 


Figure  60.  Comparison  of  Minimum,  Maximum,  and  Chassis  Potentials  as  a  Function  of 
Time  Computed  in  Three  Different  Ways  for  6  a.m.  on  January  1,  2000 
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Midnight:  January  1,2000;  NASA  Worst  Case  environment 


—m—  Handbook -Chassis 
— Hanoamk-M  in 
Han^feook-Max 
Maxwell  ian-Chassis 
Maxwell  atvMin 
— ■ —  Maxwell  an-Max 
—l —  Table-Chassis 

- Table-M  n 

- Table-Max 


Figure  61.  Comparison  of  Minimum,  Maximum,  and  Chassis  Potentials  as  a  Function  of 
Time  Computed  in  Three  Different  Ways  for  Midnight  on  January  1, 2000 

SAM:  January  1,  2000;  NASA  Worst  Case  environment 


■  Handbook -Chassis 
♦  Handccok-M  in 
Handbook -M  ax 
Maxwell  -arvCfoassis 
— n—  Maxwell  a^Min 
— »—  Maxwell  tan-Max 
— i — -  Table-Chassis 

- Table-M-n 

- Table-M  aix 


Figure  62.  Comparison  of  Minimum,  Maximum,  and  Chassis  Potentials  as  a  Function  of 
Time  Computed  in  Three  Different  Ways  for  6  a.m.  on  January  1,  2000 
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8.5  Details  of  Selected  Classes 


Nascap-2k  RealTime  Prototype  is  pure  java  and  is  based  on  the  code  used  in  the  SEE  Interactive 
Spacecraft  Charging  Handbook.  The  code  that  performs  the  calculations  is  contained  in  the 
packages  BEM,  lapack  and  SEE.  The  main  driver  code  that  handles  reading  and  writing  of  files, 
setting  up  the  calculations,  etc.,  is  in  the  package  com.saic. charging.  The  code  also  uses  some 
classes  in  the  SAIC  libraries  Utils,  SpaceXML,  and  MxOrbits.  The  entire  source  is  included,  so 
the  code  can  be  modified  and  new  jars  can  be  built. 

The  main  driver  class  is  com.saic.charging.ChargeMain.  The  method  run()  handles  initialization 
and  time  stepping  of  the  calculations. 

The  configuration  file  is  read  in  and  its  data  is  stored  in  the  class  SatelliteConfigFileReader. 

The  environments  are  read  using  the  EnvironmentMapReader  class.  The  method  that  handles 
reading  in  the  tabular  data  is  EnvironmentMapReader. getEnvironmentMapFromTextFile().  This 
method  reads  data  from  an  MSM  file  and  produces  a  map  that  contains  EnvironmentData  objects 
calculated  from  the  data  on  the  file  with  the  dates  of  the  environments  as  the  keys  for  the  map. 

This  method  is  strongly  tied  to  the  format  of  the  data  supplied  in  the  environment  file  and  this 
method  will  need  to  be  changed  with  any  change  to  that  format. 

The  EnvironmentData  objects  can  hold  either  a  Maxwellian  or  a  tabular  environment  along  with 
the  position  and  additional  environment  data  that  is  stored  in  the  java  class 
ExtraEnvironmentData.  The  ExtraEnvironmentData  class  currently  has  the  Eclipse  indicator  that 
is  used  in  the  charging  calculation  along  with  the  lsh  and  bbO  values  that  are  read  from  the  MSM 
files,  but  not  used  in  the  calculations.  The  EnvironmentData  class  has  a  method 
hasSameValuesAs()  that  is  used  to  determine  if  two  environments  are  the  same.  This  method 
may  need  to  be  modified  if  there  are  changes  to  the  format  of  the  Environment  data.  It  also  has  a 
method  writeToXml()  to  write  the  environment  to  the  output  data  files. 

The  class  SystemStateRecorder  handles  writing  the  state  of  the  system  to  the  xml  files  and  can  be 
changed  to  modify  what  is  written. 

The  SUR  file  contains  the  data  that  the  real-time  charging  code  provides  to  the  analysis  system. 

It  is  written  in  the  method  writeSURFile()  of  the  main  class  ChargeMain.  The  file  name,  data 
written  and  format  of  the  data  are  as  described  in  Interface  Control  Document  for  The  SEEFS 
Satellite  Charging/Discharging  Product  dated  May  11,  2005.  The  data  in  the  SUR  file  that  is 
calculated  by  the  real-time  charging  code  are  the  chassis  potential,  the  largest  magnitude  positive 
and  negative  potential  differences  from  the  chassis  potential  and  the  confidence  levels. 


9.  MEO  RADIATION 

We  performed  an  extensive  analysis  of  the  data  from  the  CRRES  MEA  and  HEEF  electron 
detectors,  obtaining  pitch-angle  distributions  of  the  form  j(a)  =  f90  sin”  a  ,  where: 
a  =  particle  pitch  angle,  degrees 
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j(a)  =  directional  flux  at  pitch  angle  a  (e.g. ,  in  #/cm  -s-sr-keY) 
j9o  =  directional  flux  at  a  pitch  angle  of  90°;  also  called  jperp. 
n  =  anisotropy  factor 

The  parameters  /90  and  n  were  obtained  by  performing  a  least-squares  fit  to  the  measured  pitch- 
angle  distributions.  In  performing  these  fits,  we  also  identified  cases  where  (a)  the  measured 
distribution  was  classified  as  a  “butterfly”  distribution,  i.e.,  the  flux  measured  near  a  ~  45°  was 
higher  than  that  near  90°;  and  (b)  the  measured  distribution  could  not  be  adequately  described  as 
a  butterfly  distribution  or  a  sin”  a  distribution. 

These  pitch-angle  fits  were  performed  on  one-minute  averages  of  the  raw  data.  Once  the  one- 
minute  averages  were  fit,  the  fits  were  binned  in  L.  In  order  to  investigate  the  effect  of 
geomagnetic  conditions,  the  fits  were  also  binned  in  Kp.  The  results  of  these  statistically 
averaged  and  binned  fits  are  shown  in  Figure  63  to  Figure  66.  These  figures  show  the  HEEF 
measurements  at  1.6  MeV  and  the  MEA  measurements  at  1.58  MeV.  The  following  observations 
are  made  from  these  figures: 

•  For  both  HEEF  and  MEA,  n  tends  to  increase  with  geomagnetic  activity 

•  For  both  instruments,  jperp  is  weakly  dependent  on  geomagnetic  activity,  peaking  for  Kp 
of  about  3-4. 

•  For  both  /pcrp  and  n,  both  instruments  agree  reasonably  well  for  L  >  4.  For  L  <  3.5,  the 
two  instruments  give  very  different  results. 

Figure  67  and  Figure  68  compare  the  HEEF  and  MEA  average  n  and jpeTp  as  a  function  of  energy 
at  one  L  value.  The  energy  ranges  of  the  two  instruments  overlap  from  approximately  0.65  to  1.6 
MeV.  In  this  overlap  range,  the  anisotropy  parameter  n  is  consistent  between  the  two 
instruments,  within  the  uncertainty  in  the  averages,  but  the  n  derived  for  HEEF  increases  with 
energy,  while  that  for  MEA  remains  fairly  constant.  The  energy  spectrum  in  Figure  68  shows 
that  the  value  of jperp  derived  from  the  two  instruments  is  consistent  within  the  range  of  overlap, 
but  that  the  flux  drops  off  quickly  with  energy  at  energies  above  the  region  of  overlap.  This  trend 
may  be  consistent  with  other  measurements  and  models. 

Figure  63  to  Figure  68  show  results  of  binned  and  averaged  PAD  fits.  Figure  69  to  Figure  72 
show  details  of  individual  1 -minute  average  PADs.  Figure  69  and  Figure  70  compare  directional 
fluxes  and  PAD  fits  for  MEA  and  HEEF  at  two  energies.  In  these  figures,  the  symbols  represent 
the  measured  directional  flux,  and  the  pink  lines  show  the  PAD  fit. 

The  fitting  parameters  are  given  at  the  top  of  the  figure.  Figure  71  shows  a  “typical”  butterfly 
PAD  (in  this  case  measured  by  MEA),  and  Figure  72  shows  a  typical  PAD  which  was  rejected 
by  the  fitting  procedure. 
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Figure  63.  Average  Anisotropy  Parameter  n  as  a  Function  of  L  and  Kp  for  HEEF  1.6-MeV 
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Figure  64.  Average  Anisotropy  Parameter  n  as  a  Function  of  L  and  Kp  for  MEA  1.58-MeV 

channel 
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Figure  65.  Average  jperp  as  a  Function  of  L  and  Kp  for  HEEF  1.6-MeV  Channel 
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Figure  66.  Average  jperp  as  a  Function  of  L  and  Kp  for  MEA  1.58-MeV  Channel 
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Figure  67.  Average  Anisotropy  Parameter  as  a  Function  of  Energy  for  HEEF  and  MEA 


Figure  68.  Average  jperp  as  a  Function  of  Energy  for  HEEF  and  MEA 
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Figure  69.  Comparison  of  “Typical”  Pitch-Angle  Distributions  for  MEA  (Left)  and  HEEF 

(Right),  ~  0.67  MeV 
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Figure  70.  Comparison  of  “Typical”  Pitch-Angle  Distributions  for  MEA  (Left)  and  HEEF 

(Right),  ~  0.96  MeV 
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71.  “Typical”  Butterfly  Distribution  (from  ME  A  0.69  MeV  Channel) 
Pitch  Angle  Distribution 
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Figure  72.  “Typical”  Rejected  Pitch-Angle  Distribution 
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10.  SELF-CONSISTENT  DSX  CALCULATIONS 


The  DSX  dipole  antenna  transmitter  has  two  elements,  each  several  centimeters  in  diameter  and 
40  meters  long.  Due  to  the  ease  of  electron  collection  by  positive  objects,  the  positive  element  is 
at  near  zero  or  small  positive  potential  relative  to  the  ambient  ionosphere,  while  nearly  the  full 
applied  potential  appears  on  the  negative  element.  Because  the  frequencies  of  interest  are 
comparable  to  the  ion  plasma  frequency,  the  sheath  structure  is  at  some  intermediate  state 
between  the  “ion  matrix”  or  “frozen  ion”  limit  (which  assumes  the  ions  are  stationary  and 
contribute  ambient  ion  density  to  the  space  charge)  and  the  equilibrium  space  charge  limit  (in 
which  the  ions  assume  a  steady-state  space  charge  limited  distribution  of  charge  and  current). 
Thus,  calculation  of  the  sheath  structure  and  of  ion  collection  by  the  antenna  elements  requires 
dynamic  (specifically,  particle-in-cell,  PIC)  treatment,  at  least  for  the  ions.  Nascap-2k  can  be 
used  to  perform  all  four  simulations  of  interest:  (1)  equilibrium  space  charge  sheath;  (2)  “frozen 
ion”  sheath;  (3)  dynamic  PIC  ions  with  fluid  (Boltzmann  or  barometric)  electrons  (hybrid  PIC); 
(4)  dynamic  PIC  ions  and  electrons  (Full  PIC).  One-  and  three-dimensional  Full  PIC  calculations 
[3]  and  three-dimensional  hybrid  PIC  calculations  [4  &  5]  were  reported  at  Spacecraft  Charging 
Technology  Conferences  in  2003,  2005  and  2007.  Recent  improvements  in  Nascap-2k  allow  us 
to  perform  higher  fidelity  simulations  and  better  explore  their  implications  than  was  previously 
possible.  The  calculations  reported  here  use  an  updated  geometric  model,  have  a  larger 
computational  space,  include  a  larger  number  of  macroparticles,  and  have  improved  resolution 
about  the  antenna. 

The  plasma  is  modeled  using  the  hybrid  PIC  approach  with  PIC  ions  and  fluid  barometric 
electron  densities.  The  plasma  response,  collected  ion  currents,  and  chassis  floating  potential  are 
computed  self-consistently  with  a  near-square-wave  bias  applied  to  the  antenna  elements. 

Particle  injection  and  splitting  are  used  to  replenish  the  plasma  depleted  at  the  boundary, 
represent  the  thermal  distribution,  and  maintain  appropriately  sized  macroparticles.  Limitation  of 
current  due  to  the  thermal  distribution  of  ions  and  the  resulting  angular  momentum  barrier  are 
included.  The  different  plasma  responses  above  and  below  the  ion  plasma  frequency  are 
discussed.  The  ion  current  density  is  also  computed.  Section  7  discusses  the  computation  of  the 
transverse  surface  current  and  the  volume  electron  current. 

10.1  Geometry  and  Grid 

The  DSX  spacecraft  shown  in  Figure  73  consists  of  an  EELV  Secondary  Payload  Adapter 
(ESP A)  ring  to  which  the  Avionics  Module  and  the  Payload  Module  are  attached.  The  Avionics 
Module  is  essentially  the  spacecraft  bus  and  includes  a  deployable  solar  array.  The  Payload 
Module  contains  most  of  the  DSX  specific  components.  Of  interest  here  are  the  Y  and  Z  antenna 
masts,  both  built  by  ATK  Space  Systems  Inc.  The  Y  antenna  is  80  meters  in  length  (tip-to-tip) 
and  functions  as  a  VLF  receive  and  transmit  antenna.  The  DSX  Z  antenna  is  16  meters  in  length 
(tip-to-tip)  and  functions  as  a  VLF  receive  antenna  in  a  cross-dipole  configuration  with  the  Y 
antenna.  The  Y  antenna  boom  is  a  truss  consisting  of  Graphite-Epoxy  (Gr/Ep)  longerons  and 
batten  elements  with  steel  diagonals.  Copper  wire  is  run  the  full  length  of  each  truss’s  three 
longerons,  attached  at  every  other  joint.  The  Z  antenna  boom  is  a  similar  truss  with  S-2  glass 
(fiberglass)  material  for  the  longerons  and  battens  instead  of  the  Gr/Ep.  [6] 
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Figure  73.  DSX  Spacecraft 

The  Nascap-2k  model  used  is  shown  in  Figure  74  and  Figure  75.  The  antenna  elements  are  six- 
sided,  0. 1  m  diameter  cylinders.  This  diameter  was  chosen  to  match  the  capacitance  of  the  ATK 
mesh  boom.  The  antenna  elements  are  76%  transparent  to  match  the  collecting  area  of  the  ATK 
booms.  The  solar  array  is  4.2  m  by  1.15  m  by  0.1  m.  The  18x34x18  volume  grid  is  shown  in 
Figure  76  and  Figure  77.  The  mesh  unit  of  the  outermost  grid  is  4.72  m,  and  the  mesh  unit  of  the 
innermost  grid  is  14.75  cm. 
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Figure  75.  Expanded  View  of  Center  Portion  of  Nascap-2k  Model  of  DSX 


Figure  76.  Grid  Used  for  DSX  Calculations 
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Figure  77.  Close-up  of  Center  of  Grid  Used  for  DSX  Calculations. 

10.2  Calculation  Parameters 

The  parameters  of  the  calculations  are  shown  in  Table  4.  For  the  low  density  case  (Case  1)  we 
simulated  the  applied  frequency  above  the  ion  plasma  frequency  only.  For  the  high  density  case 
we  simulated  the  system  both  above  (Case  2)  and  below  (Case  3)  the  ion  plasma  frequency.  One 
antenna  element  is  floating  and  the  other  has  a  variable  bias  with  respect  to  the  first.  The  other 
components  all  float  together.  Floating  potentials  adjust  automatically  to  account  for  incident 
plasma  current.  The  applied  bias  consists  of  the  first  two  Fourier  components  of  a  square  wave 
with  amplitude  of  1  keV  and  the  indicated  frequency.  The  shape  is  shown  in  Figure78.  The 
timestep  is  set  so  that  there  are  50  timesteps  per  cycle. 
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Table  4.  Parameters  of  Calculations  Shown 


Case  1 

Case  2 

Case  3 

Density  (m'3) 

108 

109 

109 

Temperature  (eV) 

1 

1 

1 

Species 

H+ 

H+ 

H+ 

Plasma  frequency  (kHz) 

2 

6.6 

6.6 

Frequency  (kHz) 

10 

12 

2 

Splitting  of  initial  macroparticles 

All  grids 

All  grids 

All  grids 

Splitting  on  subgrid  entry 

When  needed 

When  needed 

When  needed 

Macroparticle  injection 
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Figure  78.  Applied  Bias  Values  and  Resulting  Antenna  Element  Potentials  in  the  Absence 

of  Plasma 
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10.3  Particle-in-Cell  Techniques 

Nascap-2k  has  primarily  been  used  to  model  quasi-static  phenomena.  However,  there  are  many 
physical  processes  of  interest  whose  timescales  require  a  dynamic  approach  such  as  a  PIC 
technique.  Examples  of  such  processes  are  breakdown  phenomena,  plasma  kinetics,  and  sheath 
structure  about  surfaces  with  potentials  that  change  on  a  timescale  comparable  to  the  time  it 
takes  an  ion  to  cross  the  sheath.  PIC  techniques  can  also  be  used  to  address  problems  in  which 
analytic  representations  of  the  environmental  currents  are  inadequate,  such  as  in  a  spacecraft 
wake  or  in  a  cavity.  A  steady-state  PIC  technique,  in  which  the  ion  space  charge  density  is 
computed  from  macroparticles  tracked  from  the  boundary  of  the  computational  space  until  they 
are  collected  or  exit  the  computation  space,  was  successfully  used  to  model  the  CHAWS 
experiment.  [7]  In  addition,  PIC  techniques  can  be  useful  when  developing  analytic  models.  To 
facilitate  these  modeling  techniques,  the  ability  to  perform  various  types  of  PIC  calculations  was 
built  into  Nascap-2k. 

In  a  hybrid  PIC  calculation,  the  problem  is  initialized  by  creating  ion  macroparticles  throughout 
the  grid  to  represent  a  constant  particle  density.  The  ion  macroparticles  are  tracked  for  one 
timestep  and  then  volume  potentials  are  computed  using  the  resulting  ion  density  and  a 
barometric  electron  density.  The  process  is  repeated  for  the  time  period  of  interest.  In  a  full  PIC 
calculation  both  electron  and  ion  macroparticles  are  tracked  and  volume  potentials  are  computed 
using  the  resulting  plasma  density. 

The  barometric  potential  is  defined  by  4>b  =  0  ln(pj/en),  where  n  and  9  are  the  ambient  plasma 
density  and  temperature  and  pi  is  the  local  ion  density  (obtained  from  tracking).  At  plasma 
potentials  below  0b,  the  electron  density  is  the  barometric  electron  density  and  above  0b,  the 
electron  density  linearly  increases.  Charge  stabilization  is  then  applied.  [8] 
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To  replace  macroparticles  that  are  collected  by  the  probe  or  escape  from  the  grid,  it  is  necessary 
to  periodically  inject  macroparticles  from  the  boundary.  This  allows  for  the  calculation  of  current 
for  longer  time  periods.  In  hybrid  PIC  calculations  without  boundary  injection,  the  low  field 
region  near  the  boundary  of  the  problem  develops  a  significant  negative  potential  due  to  ion 
depletion.  Boundary  injection  keeps  these  potentials  near  zero  by  replenishing  the  ions  that  have 
been  collected  or  escaped. 

Macroparticles  may  be  split  when  they  enter  a  more  finely  resolved  region  or  when  they  are 
created  either  at  the  boundary  or  throughout  the  volume  at  problem  initialization.  There  are  two 
major  reasons  for  splitting  macroparticles:  one  physical  and  one  numeric.  Even  at  moderate 
potentials,  thermal  effects  can  reduce  collected  currents.  Some  particles  near  the  sheath  edge 
have  enough  thermal  velocity  perpendicular  to  the  electric  field  that  angular  momentum 
conservation  prohibits  collection.  Particle  splitting  allows  for  a  representation  of  the  thermal 
distribution  in  the  initial  particle  distribution  and  in  particles  injected  from  the  boundary.  From  a 
numeric  point  of  view,  particle  splitting  can  be  used  to  keep  the  particle  weight  appropriate  to  the 
grid  size  and  to  help  maintain  the  smoothness  of  the  distribution.  A  large  particle  that  originated 
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in  an  outer  grid  is  split  in  velocity  in  such  a  way  that  it  preserves  the  plasma  temperature  and 
eventually  becomes  distributed  over  several  volume  elements  in  an  inner  grid.  Because  high- 
field  regions  are  often  of  interest,  particles  are  split  in  velocity  space  only,  as  spatial  splitting 
would  raise  problems  with  energy  conservation.  The  details  are  discussed  in  Reference  5. 

10.4  Results 

The  time  history  of  the  surface  potentials  and  tracked  ion  currents  to  surfaces  for  the  three 
calculations  are  shown  below  in  Figure  79  through  Figure  87.  In  each  calculation,  there  is  an 
initial  transient  in  the  potentials  and  currents.  In  the  low-density  Case  1,  it  takes  about  6  cycles 
until  the  potential  variation  of  the  antenna  elements  settles  down  to  approximate  a  square  wave 
of  amplitude  500  V  about  a  value  slightly  more  negative  than  -500  V.  In  the  higher  density  cases, 
it  takes  about  one  cycle.  In  all  cases,  the  body  goes  negative  to  about  -200  V  and  then  returns  to 
zero  at  a  rate  that  depends  on  the  plasma  density. 

When  the  excitation  frequency  is  above  the  ion  plasma  frequency  (Case  1  and  Case  2)  the  plasma 
current  to  an  antenna  element  increases  sharply  when  a  potential  is  applied  to  the  element,  and 
continues  substantially  into  the  unbiased  half-cycle,  as  seen  in  Figure  81  and  Figure  84.  Below 
the  ion  plasma  frequency  (Case  3)  the  ion  current  peaks  sharply  when  the  potential  is  applied, 
falls  off  substantially  during  the  half-cycle,  and  drops  sharply  to  zero  when  the  potential  is 
removed,  as  seen  in  Figure87.  The  difference  can  be  understood  in  terms  of  the  sheath  being 
depleted  of  ions  within  every  cycle  at  low  frequency,  but  not  at  high  frequency.  For  excitation 
frequencies  below  the  plasma  frequency,  the  initial  burst  of  ion  current  depletes  the  sheath.  The 
lower  current  levels  during  the  second  half  of  the  biasing  period  are  a  result  of  ions  being 
attracted  from  the  sheath  edge.  This  shows  up  in  the  phase  shift,  tabulated  in  Table  4.  The 
negative  phase  shift  for  Case  3  is  a  result  of  the  current  maximum  occurring  shortly  after  the 
onset  of  the  high  voltage,  while  the  positive  phase  shifts  for  cases  1  and  2  show  the  delay  in 
current  response  to  the  high  voltage  application.  These  phase  shifts  are  much  less  than  the  ninety 
degree  phase  shift  expected  in  the  linear  regime. 
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Figure  79.  Time  Dependence  of  Antenna  Element  Potentials  for  Case  1 
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Figure  80.  Potential  and  Collected  Ion  Current  of  Antenna  Element  2  for  Case  1 
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Figure  81.  Potential  and  Collected  Ion  Current  of  Antenna  Element  2  for  Case  1 

(Expanded  Scale) 
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Figure  82.  Time  Dependence  of  Antenna  Element  Potentials  for  Case  2 
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Figure  83.  Potential  and  Collected  Ion  Current  of  Antenna  Element  2  for  Case  2 
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Figure  84.  Potential  and  Collection  Ion  Current  of  Antenna  Element  2  for  Case  2 

(Expanded  Scale) 
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Figure  85.  Time  Dependence  of  Antenna  Element  Potentials  for  Case  3 
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Figure  86.  Potential  and  Collected  Ion  Current  of  Antenna  Element  2  for  Case  3 
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Figure  87.  Potential  and  Collected  Ion  Current  of  Antenna  Element  2  for  Case  3 

(Expanded  Scale) 

Figure  88  shows  the  sheath  structure  in  the  low  density  plasma  when  the  potential  on  one 
antenna  element  is  at  its  negative  peak  at  -1.1  kV  and  the  other  is  positive  at  +26  V.  At  this 
potential  and  density,  the  sheath  is  nearly  cylindrical,  with  a  spherical  end  cap.  The  sheath 
envelops  the  positive  antenna  element,  providing  a  barrier  to  electrons,  and  thus  allowing  the 
system  to  float  more  positive  than  might  be  expected.  Presently,  this  barrier  effect  on  the  current 
collection  is  not  modeled.  In  the  immediate  vicinity  of  the  antenna,  the  density  is  elevated  as  the 
ions  are  attracted  to  the  negative  potential.  This  convergence  leaves  a  slight  decrease  in  density 
at  the  edge  of  the  sheath.  The  elevated  density  around  the  positive  element  is  left  over  from  the 
previous  half  cycle. 

Figure  89  shows  the  sheath  structure  when  both  antenna  elements  are  at  the  same  potential,  near 
-570  V.  The  sheath  is  slightly  larger  about  the  top  element,  which  had  previously  been  near 
-1  kV,  while  a  large  volume  of  the  lower,  previously  positive,  sheath  remains  near  ambient  ion 
density  due  to  ion  inertia  effects.  In  the  immediate  vicinity  of  the  top  antenna  element,  the  ion 
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density  is  enhanced  by  an  order  of  magnitude  due  to  ions  that  orbit  the  probe  but  are  not 
collected.  However,  throughout  most  of  the  sheath  the  ion  density  is  depressed,  leading  to  less 
shielding  and  a  larger  sheath. 


Figure  88.  Sheath  Structure  (Potentials  and  Densities)  for  Low  Density  Calculation  at 
2.062  ms  -  Antenna  Elements  are  at  -1173  V  And  +26  V,  the  Largest  Differential  of  the 

Cycle  -  Ambient  Plasma  Density  is  1.8  V  m'2 
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Figure  89.  Sheath  Structure  (Potentials  and  Densities)  for  Low  Density  Calculation  at 
2.05  ms  -  Antenna  Elements  are  at  -575  V  and  -573  V  -  the  Bottom  Antenna  Element  was 
Previously  Near  Plasma  Ground,  while  the  Top  Element  was  Recently/Previously  Near  -1 

kV  -  Ambient  Plasma  Density  is  1.8  V  m'2 
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The  electromagnetic  waves  from  the  antenna  are  created  by  current  flow.  The  near  field  current 
has  three  components:  transverse  surface  current  along  the  antenna  elements,  plasma  ions,  and 
plasma  electrons.  The  plasma  ion  currents  are  obtained  directly  from  the  PIC  calculations,  and 
are  discussed  below.  The  surface  currents  and  volume  electron  currents  are  obtained  via  pseudo¬ 
potential  methods  as  discussed  in  Section  7. 

Figure  90,  Figure  91,  and  Figure  92  show  the  X  and  Y  components  of  the  volume  ion  current 
within  the  computational  grid  when  the  antenna  elements  are  at  the  same  potential  and  at  one- 
fifth  and  three-fifths  of  a  half-cycle  later  (near  the  first  and  second  peaks  of  the  applied 
potential).  The  left  panels  of  these  figures  illustrate  how  ions  that  are  accelerated  toward  the 
antenna  and  orbit  around  it  form  periodic  outwardly  moving  blocks  of  ions.  The  details  are  as 
follows: 

(1)  Adjacent  to  the  upper  antenna  element  in  Figure  90  (which  is  in  the  process  of  switching 
from  negative  to  positive)  we  see  the  ions  (which  had  been  orbiting  the  antenna)  begin  to 
move  outward  as  the  attractive  force  is  removed. 

(2)  In  Figure  91,  with  the  electric  field  now  slightly  repulsive,  a  solid  block  of  outwardly 
moving  ions  is  seen  next  to  the  upper  antenna  element. 

(3)  In  Figure  92  the  block  of  outwardly  moving  ions  separates  from  the  upper  antenna 
element. 

(4)  Returning  to  Figure  90  and  now  focusing  on  the  lower  antenna  element  (which  is  a  half¬ 
cycle  removed  from  the  upper)  we  see  that  the  outwardly  moving  block  has  radially 
expanded.  The  inner  portion  of  the  block  consists  of  relatively  slowly  moving  ions  that 
will  be  recaptured  by  the  increasing  negative  potential.  The  outer  portion  of  the  block, 
which  consists  of  more  energetic  ions  that  have  escaped  to  the  outer,  weak  field  region  of 
the  sheath,  will  continue  to  move  outward.  In  the  high  field  region  immediately  adjacent 
to  the  newly  negative  antenna  element,  the  current  is  high  and  ion  motion  consistent  with 
the  field  direction. 

(5)  Figure  91  and  Figure  92  show  continued  outward  motion  of  the  block  around  the  lower 
antenna  element.  Closer  to  the  antenna,  ions  are  moving  inward  in  preparation  for 
forming  the  next  block. 

(6)  Returning  to  the  upper  antenna  element,  the  block  we  have  been  following  appears  (a  full 
cycle  later)  in  Figure  90  about  halfway  between  the  antenna  and  the  problem  boundary. 

Its  outward  motion  continues  in  Figure  91  and  Figure  92. 

The  right  panels  of  these  figures  show  that  the  Y-component  of  ion  current  near  the  outboard 
three-quarters  of  each  element  is  dominated  by  ions  flowing  inward  along  the  antenna  after  being 
attracted  by  the  field  near  the  boom  tip.  The  Y-component  of  ion  current  for  roughly  the  inboard 
quarter  of  each  element  is  flowing  upward  in  Figure  90,  toward  the  previously  negative  element. 
Due  to  ion  inertia,  the  current  near  the  root  of  the  upper  antenna  element  is  still  flowing  upward 
in  Figure  91,  even  though  the  antenna  polarity  has  reversed.  After  twenty  milliseconds  of 
applying  the  downward  field  (Figure  92)  the  ion  flow  has  reversed  to  be  consistent  with  the 
applied  field,  and  is  now  downward  near  the  antenna  element  roots. 
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Figure  90.  X  and  Y  Components  of  Volume  Current  Density  Due  to  Ions  at  2.55  ms  - 
Antenna  Elements  are  at  -576  V  and  -578  V  -  Previously,  the  Bottom  Antenna  Element 
was  Near  Plasma  Ground  and  the  Top  Strongly  Negative 
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Figure  91.  X  and  Y  Components  of  Volume  Current  Density  Due  to  Ions  at  2.562  ms  - 
Antenna  Elements  are  at  -1152  V  (Lower)  and  1.8  V  (Upper) 
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Figure  92.  X  and  Y  Components  of  Volume  Current  Density  Due  to  Ions  at  2.582  ms  - 
Antenna  Elements  are  at  -1058  V  (Lower)  and  0.3  V  (Upper) 

Figure  93  and  Figure  94  show  the  vector  potential,  rate  of  change  in  the  vector  potential,  and  the 
magnetic  field  resulting  from  the  transverse  surface  current.  The  radiating  electric  field  is 
proportional  to  the  rate  of  change  of  the  vector  potential. 
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Figure  93.  Y  Component  of  the  Vector  Potential  and  Rate  of  Change  of  the  Vector 
Potential  at  2.05  ms  as  a  Result  of  the  Transverse  Currents  Shown  in  Figure  50 
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Figure  94.  Z  Component  of  the  Magnetic  Field  at  2.05  ms  as  a  Result  of  the  Transverse 

Currents  Shown  in  Figure  50 
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10.5  Magnetic  Field 

We  also  looked  at  the  effect  of  the  magnetic  field  on  the  ion  dynamics  calculations  by  repeating 
the  10  kHz  calculation  with  a  6  pT  (0.06  gauss)  field  in  the  Y  (parallel)  direction  and  in  the  X 
(perpendicular)  direction.  For  this  magnetic  field  and  1  kV,  the  gyroradius  of  hydrogen  ions  is 
760  m.  In  both  cases,  we  started  the  calculation  with  the  result  with  no  magnetic  field  at  1  ms.  As 
expected  there  are  no  discemable  differences  in  the  ion  surface  and  volume  currents  with  the 
magnetic  field  included  in  the  tracking  calculations. 

10.6  5  kV 

We  started  a  10  kHz,  5  kV  calculation.  The  time  histories  of  the  potentials  and  currents  are 
shown  in  Figure  95  through  Figure  97.  Other  than  a  factor  of  5  increase  in  the  voltage  and 
collected  current,  at  1.3  ms,  the  results  are  similar  to  those  at  1  kV  at  1.3  ms.  The  sheath  structure 
near  1.3  ms  for  1  kV  and  5  kV  is  shown  in  Figure  98  and  Figure  99.  Other  than  the  much  larger 
sheath  at  5  kY,  the  sheath  is  much  the  same.  The  ion  currents  at  1.3  ms  are  shown  in  the  final 
figure.  At  the  higher  potential,  the  released  ion  current  may  be  more  concentrated  toward  the 
center  of  the  antenna. 
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Figure  95.  Time  Dependence  of  Antenna  Element  Potentials  for  5  kV  at  10  kHz 
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Figure  96.  Potential  and  Collection  Ion  Current  of  Antenna  Element  2  for  5  kV  at  10  kHz 
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Figure  97.  Potential  and  Collection  Ion  Current  of  Antenna  Element  2  for  5  kV  at  10  kHz 
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Figure  98.  Sheath  Structure  (Potentials)  at  1.298  ms  for  5  kV  (Left)  at  10  kHz  and  1  kV 

(Right)  at  10  kHz 


Figure  99.  Sheath  Structure  (Density)  at  1.298  ms  for  5  kV  (Left)  at  10  kHz  and  1  kV 
(Right)  at  10  kHz  -  Ambient  Plasma  Density  is  1.8  V  m'2 
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10.7  Conclusions 


We  have  used  Nascap-2k  to  study  the  plasma  interactions  of  a  high  voltage  near-square-wave 
VLF  antenna  in  MEO  plasma.  The  geometric  model  captures  the  main  features  of  the  DSX 
spacecraft.  The  antenna  elements  are  biased  ±1000  V,  at  frequencies  above  or  below  the  ion 
plasma  frequency.  For  these  very  high  applied  voltages,  the  plasma  response  is  nonlinear,  and 
the  sheaths  are  as  much  spherical  as  cylindrical,  especially  at  lower  densities.  In  modeling  the 
system  it  is  important  to  use  particle  splitting  and  injection  techniques  that  replenish  depleted 
plasma  at  boundaries,  maintain  appropriately  sized  macroparticles,  and  provide  a  reasonable 
representation  of  the  plasma  thermal  distribution.  When  this  is  done,  the  incident  current  is  in 
reasonable  agreement  with  orbit-limited  predictions. 

The  sheath  is  depleted  of  ions  within  every  cycle  when  the  antenna  is  driven  below  the  ion 
plasma  frequency,  but  not  at  high  frequency.  The  current  lags  the  applied  voltage  for  excitation 
frequencies  above  the  ion  plasma  frequency  by  1 1°  at  the  low  density  and  13°  at  the  high 
density,  which  is  much  less  than  the  ninety  degree  shift  expected  in  the  linear  regime.  At  low 
frequency  (calculated  only  for  the  high  density  case)  the  current  leads  the  applied  voltage  by  13°. 

There  are  three  near-field  sources  of  current:  current  that  flows  along  the  antenna  element 
surfaces  (supplied  by  the  power  supply),  plasma  ion  currents  that  are  attracted  and  released  by 
the  applied  potential,  and  currents  of  plasma  electrons  that  must  vacate  an  expanding  sheath  and 
refill  a  contracting  one.  The  plasma  ion  current  is  computed  during  the  PIC  plasma  calculation, 
and  has  been  discussed  here.  We  predict  the  periodic  launch  of  blocks  of  energetic  ions  radially 
outward  into  the  cold  plasma.  We  have  also  developed  a  pseudo-potential  method  for  calculating 
the  surface  currents  and  plasma  electron  currents  self-consistently  with  the  Nascap-2k  PIC 
calculation  of  the  sheath  dynamics;  those  results  are  described  in  Section  7.  These  methods  allow 
us  to  quantitatively  predict  the  disturbances  to  the  presheath  plasma  that  result  from  the 
dynamics  of  the  highly  nonlinear  sheath  region,  enabling  the  use  of  linear  plasma  theory  to 
predict  generation  of  propagating  plasma  excitations. 

Improvements  to  Nascap-2k  over  the  past  six  years  provide  the  capabilities  needed  to  predict  the 
near  field  behavior  of  electromagnetic  transmission  from  a  VLF  transmission  in  a  plasma.  The 
first  of  these  improvements  is  the  enhancement  of  Nascap-2k  PIC  capabilities  (particle  splitting 
and  injection)  and  integration  of  the  PIC  algorithms  with  earlier  steady-state  algorithms.  The 
replacement  of  the  original  (1980s)  database  with  a  more  flexible  one  with  orders  of  magnitude 
more  storage  capability  removes  the  limitation  on  the  size  and  complexity  of  problems  that  can 
be  addressed  and  also  facilitates  the  use  of  algorithms  that  require  the  storage  of  more 
information  about  the  local  plasma. 


11.  CONCLUSION 

The  program  successfully  improved  the  plasma  engineering  capability  available  to  the  spacecraft 
community.  New  algorithms  for  Nascap-2k  were  developed,  incorporated,  tested,  and  validated. 
Nascap-2k' s  PIC  computational  capabilities  were  dramatically  improved.  Nascap-2k' s  charging 
capabilities  were  enhanced.  A  transverse  surface  current  model  was  added  and  an  algorithm 
developed  to  compute  the  volume  electron  current.  Improvements  were  also  made  to  the 
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usability  of  the  software.  The  new  database  and  memory  manager,  N2kDB,  enables  further 
Nascap-2k  development  and  allows  larger  calculations  to  be  performed. 

The  prototype  of  Nascap-2k  RealTime  is  being  evaluated  by  the  SEEFS  program. 

A  set  of  self-consistent  calculations  of  the  plasma  response  to  a  high  voltage  square  wave  VLF 
DSX  antenna  were  performed  in  2007  and  repeated  with  higher  fidelity  at  the  end  of  the  contract. 
These  calculations  used  the  new  and  enhanced  capabilities  developed  under  this  program. 

Much  of  the  work  performed  under  this  program  has  been  presented  at  conferences  and 
published  in  refereed  journals.  Details  of  the  work  are  documented  in  this  and  other  AFRL 
reports.  At  this  writing,  the  next  code  version,  including  the  new  database  and  memory  manager 
as  well  as  other  enhancements,  is  undergoing  beta  testing. 


The  Nascap-2k  enhancements  and  new  algorithms  enable  a  wide  variety  of  PIC  calculations  to 
be  performed.  In  addition,  the  improvements  form  the  groundwork  for  further  algorithm  and 
software  development. 
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